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Abstract 


A mixed method of approximation based on Reissner’s variational principle is 
developed for the linear analysis of interlaminar stresses in laminated composites, with 
special interest in laminates that contain terminated internal plies (dropped-ply 
laminates). Two models are derived, one for problems of generalized plane deformation 
and the other for the axisymmetric response of shells of revolution. A layerwise approach 
is taken in which the stress field is assumed with an explicit dependence on the thickness 
coordinate in each layer. The dependence of the stress field on the thickness coordinate is 
determined such that the three-dimensional equilibrium equations are satisfied by the 
approximation. The solution domain is reduced to one dimension by integration through 
the thickness. Continuity of tractions and displacements between layers is imposed. 

The governing two-point boundary value problem is composed of a system of both 
differential and algebraic equations (DAEs) and their associated boundary conditions. 
Careful evaluation of the system of DAEs was required to arrive at a form that allowed 
application of a one-step finite difference approximation. A two-stage Gauss implicit 
Runge-Kutta finite difference scheme was used for the solution because of its relatively 
high degree of accuracy. Patch tests of the two models revealed problems with solution 
accuracy for the axisymmetric model of a cylindrical shell loaded by internal pressure. 

Parametric studies of dropped-ply laminate characteristics and their influence on the 
interlaminar stresses were performed using the generalized plane deformation model. 
Eccentricity of the middle surface of the laminate through the ply drop-off was found to 
have a minimal effect on the interlaminar stresses under longitudinal compression, 
transverse tension, and in-plane shear. A second study found the stiffness change across 
the ply termination to have a much greater influence on the interlaminar stresses. 
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Correlations between the stiffness ratio of the thick to the thin sections of the laminates 
and the magnitude of a parameter based on a quadratic delamination criterion were found 
to be surprisingly good for longitudinal compression and in-plane shear loadings. For 
laminates with very stiff terminated plies loaded in longitudinal compression, inclusion of 
a short insert of softer composite material at the end of the dropped plies was found to 
significantly reduce the interlaminar stresses produced. 
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Chapter 1: Introduction 


The application of advanced composites in aerospace structures has increased 
dramatically in the past decade. Their advantages in stiffness and strength to weight ratios 
over conventional metals allow lighter, more efficient aircraft. In addition, their anisotropic 
behavior allows the design of components with response characteristics that are not 
possible with isotropic materials. The most visible example of such an application is the 
forward swept wing of the X-29 technology demonstrator which has bending-twist 
coupling characteristics that prevent the aeroelastic divergence that had previously made 
such a design impractical. 

When composite materials are used in aerospace structural components, it is usually 
desirable to tailor the material to match the localized strength and stiffness requirements in 
order to minimize the weight. For a fibrous composite laminate composed of unidirectional 
layers, this is often done by changing the number of plies. This abrupt change in thickness, 
which is referred to as a ply drop-off (other sources have used the terms tapered laminates 
and laminates with terminated plies), introduces a stress concentration which can lead to 
premature delamination failure of a laminate. It is this interlaminar stress concentration that 
is the subject of the present research. 

A common application of the ply drop-off is a composite aircraft wing skin, or 
stabilizer skin, which requires much greater thickness at the root as compared to the tip. 
When fasteners are used to attach two components together, laminates must be built up to 
handle the concentrated bearing loads. The multi-segment composite solid rocket booster 
casing for the space shuttle currendy under development is a highly visible example of this 
application. 1 - 2 Other applications of composites with dropped plies which have been 
mentioned in the literature include helicopter yokes, 3 helicopter rotor flexbeams, 4 - 5 and 
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flywheels . 6 


1.1 Literature Review 

Most of the literature on the subject of ply drop-offs has been published in the past 
decade. Research has been done on the response, failure initiation, damage tolerance, 
failure modes and even beneficial effects on failure loads of ply termination in laminated 
composites. Longitudinal and transverse as well as internal and external ply drop-offs have 
been examined under various loading conditions and environmental influences. This 
literature review will concentrate on the static response of dropped-ply laminates without 
stability considerations and will therefore exclude studies dealing with fatigue or buckling. 

1.1.1 Classification and Overview of Parametric Effects 

There are several types of ply drop-offs which differ depending upon which plies 
through the thickness are dropped and the orientation of the termination with respect to the 
loading. An external ply drop-off is one in which the dropped plies are on a surface of the 
laminate while an internal ply drop, as its name implies, has plies terminated from the 
interior. Because of the large peeling stresses present in such a design, external ply drop- 
offs are relatively weak and are therefore avoided when possible. The most common 
instance where an external ply drop-off is used is for attachment of stiffeners . 7 Internal ply 
drop-offs, which require the remaining plies to contour around the ends of the dropped 
plies, while not without their own problems, are more common in stiffness tailoring. 

A ply drop-off can be classified further depending upon the orientation of the load path 
with respect to the direction normal to the terminated edge. Excluding shear loadings, 
which have received little attention in the literature, a longitudinal ply drop-off is oriented 
such that the normal to the terminated edge is perpendicular to the load and a transverse 
drop has the normal parallel to the loading direction. 

A ply drop-off can also be classified based on symmetry. Symmetric ply drop-offs are 
those where all plies that are dropped on one side of the middle surface have a counterpart 
of the same orientation dropped symmetrically with respect to this middle surface. In 
addition, the middle surface does not shift in the thickness direction through the ply drop- 
off region. A configuration not classified as such is asymmetric and couples extension and 
bending. Aircraft wing skins typically utilize asymmetric ply drop-offs in order to have the 
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surface exposed to the airflow as flat as possible. 

External ply drops were examined by Daoust and Hoa, 3 Wu and Webber, 8 Wu, 9 
and by Kubr. 7 Wu and Webber performed a finite element analysis of external ply drops 
under transverse tensile loading and found very high peak stresses in the comer of the 
drop. When a resin fillet was added to the analysis, the peak interlaminar stresses were 
reduced to about half their original values. The paper by Wu was a follow up of this 
analysis, in which nonlinear material behavior was included to account for the 
redistribution of stresses in the resin that would occur in the presence of the peak stresses. 
Without a fillet, resin nonlinearity also had the effect of reducing the peak interlaminar 
stresses by about half. No comparison was made to a uniform thickness laminate so the 
degree of strength degradation was not determined. Daoust and Hoa used finite elements to 
examine the effects of several parameters on the strength of transverse ply drop-offs 
including external versus internal configurations. Under tension, bending and torsion they 
found that internal ply drop-offs are typically twice as strong as external ply drop-offs. 
Kubr looked at external ply drop-offs in the context of co-cured stringer reinforced 
composite plates. He used a finite element analysis and examined the stress singularity at 
the reentrant comer using a highly refined mesh in this area. 

A few papers examined longitudinal ply drop-offs. In each case, the purpose of 
introducing longitudinal drop-offs was to alter the stress field in the free edge region in an 
attempt to prevent delamination. Chan and Ochoa 10 used finite elements in analyzing 
laminates with 90° plies dropped symmetrically just inside of the free edge and found that 
under tensile loading, interlaminar normal stresses were reduced by 85% and interlaminar 
shear by 43%. In an accompanying experimental program, the tensile strength was found 
to be increased by 40% and premature delamination was eliminated. In a follow up 
study, 11 the authors also examined the longitudinal ply drop-off under bending and torsion 
loads. They found that under these loadings, the maximum interlaminar stresses were not 
altered significantly and in most cases were increased a small amount. Pogue and Vizzini 12 
looked at four methods of edge alteration in tension specimens, including dropping plies, in 
a purely experimental study. Vizzini 13 later added a finite element analysis which agreed 
reasonably well with experiments in predicting the strain level at failure initiation. These 
studies reached a different conclusion than Chan and Ochoa regarding the merits of 
longitudinal ply drops. Working with symmetric delamination-critical laminates composed 
of ±15° and 0° layers, Vizzini found that dropping plies from these laminates usually 
caused failure to initiate internally at the end of the dropped layer. The delamination 
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initiation load was found to be lowered by as much as 10% or increased by as much as 
16% relative to an unaltered specimen, depending on the configuration. This difference in 
conclusions between these two studies may be attributed to the fact that for the laminates 
studied by Chan and Ochoa, the mismatch in Poisson’s ratio was critical while in those that 
Vizzini looked at, the mismatch in the coefficients of mutual influence was critical. 

The remainder of this overview concentrates on transverse internal ply drop-offs 
as their application is more prevalent and of most interest to this study. Of these, a small 
group of sources 14 - 15 - 16 examined both geometrically symmetric and asymmetric ply drop- 
offs, allowing a comparison between the two configurations. Kemp and Johnson 14 
performed a finite element analysis of two layups of T300/5208 and found that in tension, 
geometrically symmetric ply drops have a 8-21% higher strain to initial failure than 
asymmetric drops, depending on the layup and number of plies dropped. In compression 
they found the reverse, asymmetric ply drops stronger than symmetric by 7-26%. In 
addition to the change in failure strain, they also found that the location of the predicted 
failure initiation sometimes moved due to a change in symmetry conditions. Lagace and 
Cannon 1516 used an analysis employing classical laminated plate theory in approximating 
the in-plane stress concentration of an idealized asymmetric step change in thickness, 
performed experiments, and concluded that the effect on fracture of bending induced by an 
asymmetric ply drop was insignificant. However, it is worth noting here that only tension 
was considered and the asymmetric laminate tested had only one ply out of six terminated. 
In addition, a comparison to a symmetric version of the same layup was not made. 

In experimental and finite element analysis of a symmetric ply drop, Wisnom 17,18 
investigated three factors influencing failure when internal plies are terminated in 
unidirectional glass-epoxy specimens subject to either tension or compression. 
(Unidirectional specimens were used to avoid failure by free edge delamination that is 
prevalent with many angle ply specimens.) The first factor was the stress concentration 
associated with terminating internal plies, the second was the interlaminar shear and normal 
stresses induced by the tapered geometry, and the third was the effect of the curved fibers 
adjacent to the ply drop. Straight specimens with cut central plies were subject to tension to 
assess the stress concentration associated with ply termination (no taper or curved plies). 
Straight specimens were also severely tapered, or waisted, through the thickness by 
machining with a diamond wheel (no internal ply terminations or curved plies) and 
subjected to tension. Finally, the effect of fiber curvature was investigated by means of 
inserting two plies of nylon film in the middle of eight ply tensile specimens, or two layers 


4 


Introduction 



of film adhesive in the middle of thirty-two ply compression specimens. The straight 
specimens with cut central plies failed by delamination between the discontinuous and 
continuous plies, the waisted specimens delaminated due to a stress concentration 
dominated by interlaminar shear stress, and specimens with inclusions failed due to fiber 
breakage in the grips. It was concluded that the critical parameter affecting delamination in 
tapered unidirectional composites is the strain energy release rate due to the terminating 
internal plies. 

Papers by Botting, Vizzini, and Lee 19 and Fish and Vizzini 5 are similar in that they 
examined the effect of altering the sequence in which plies are dropped. Both papers 
analyzed glass/epoxy laminates with plies dropped symmetrically in three steps from 28 to 
16 plies with a taper angle of 5.7°. The Fish and Vizzini paper was purely experimental and 
all of the plies were oriented at 0° in order to minimize free edge effects. In a multi-ply 
staggered drop, they altered two characteristics of the taper region. First, the sequence that 
the plies through the thickness were dropped was varied, giving either a staircase or 
overlapped appearance. Second, the particular plies through the thickness that were 
dropped was varied, or more specifically whether the dropped plies were adjacent to each 
other or not. The result was a four way comparison between laminates that were staircased 
or overlapped and grouped or dispersed. Bending stiffness was measured at the start and as 
damage progressed. They found that delamination onset loads varied by 38% and bending 
stiffness retention varied by 56% between configurations with a tendency for tradeoff 
between the two. The overlapped-dispersed configuration was found to provide the best 
overall performance. The work by Botting, Vizzini and Lee included a three dimensional 
finite element analysis as well as experiments in examining the effects of the longitudinal 
sequence of ply termination with all dropped plies interspersed between continuous 
sublaminates. The continuous sublaminates were either [±45^0,*]* or [ 04 /± 45 2 L and the 
dropped plies were ±45 pairs. For the [0V+45J, laminates, the best configuration was one 
in which the second dropped ply from the mid-surface was overlapped over the first at the 
thin end of the drop region. This provided a 16% strength increase over the simple staircase 
configuration. Through visual inspection and x-ray photography, the [±45 2 /0 4 ] i laminates 
were found to delaminate due to free edge effects and therefore no conclusions could be 
made regarding the effect of the configuration of the drop. 

In an analysis incorporating finite elements as well as experiments, Curry et al. 20>21 
found that strength reduction due to asymmetric ply drop-offs was greater in compression 
than in tension and that the strength reduction was proportional to the stiffness change 
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across the drop. This finding was corroborated by Trethewey et al . 22 In their finite element 
study, Daoust and Hoa 3 found that internal ply drop-offs lose very little strength in torsion. 
They also found that increasing the distance between the end of the terminated plies and the 
joining of the continuous plies, thereby reducing the taper angle, lowers the magnitude of 
the interlaminar stresses generated. Kemp and Johnson 14 corroborated this finding by 
predicting increased failure strains under those circumstances in their finite element based 
analysis. 

1.1.2 Approaches to Modeling and Analysis 

The methods of analysis employed in the literature varied depending upon the 
objectives of the studies and the resources available. Because of their relative ease of 
application to geometrically complex problems, by far the most prevalent method was finite 
elements. In fact the only studies found that did not employ finite elements were the 
previously mentioned work by Lagace and Cannon , 15 - 16 Trethewey, Gillespie and 
Wilkins , 22 and a few which were purely experimental . 1 ’ 5,23,12 Lagace and Cannon 
determined the bending deflections of dropped ply laminates modeled essentially as two 
asymmetrically end- butted beams and then determined the in-plane stresses from the stress- 
strain relations based on classical laminated plate theory. The work done by Trethewey et 
al. divided a dropped-ply laminate into six subregions and employed shear deformation 
plate theory to these subregions. A damage tolerance analysis was then performed using 
linear elastic fracture mechanics. The choice of shear deformation plate theory over 
continuum finite elements was made in order to streamline the analysis for use as a 
preliminary design tool. 

Three-Dimensional Finite Element Analyses 

Some of the studies employed fully three-dimensional finite element analyses. One of 
the earliest was by Adams, Ramkumar and Walrath 24 and included the effects of nonlinear 
material response and thermal residual stresses. However, although the model was three- 
dimensional, the lateral surfaces were constrained to remain planar, essentially imposing 
generalized plane strain conditions and not accounting for free-edge effects. In addition, 
each layer was modeled by one element through the thickness except for the dropped layer 
which had two. Also, no longitudinal mesh refinement was made at the actual ply drop 
location. Hoa, Daoust, Du and Vu-Khanh 25 refined their three-dimensional mesh at the ply 
drop by a submodeling (or zooming) technique. This approach involved successive 
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reduction and refinement of the mesh in the region of interest while retaining the results of 
the previous iteration as boundary conditions for the reduced mesh. The purpose of this 
method was to have a refined mesh in the region of large stress gradients while keeping the 
number of degrees of freedom of the solution required for each pass of the finite element 
solver within the capability of the computer that was available. In an extension of this 
work, Daoust and Hoa 3 once again used three-dimensional finite elements; however rather 
than using submodeling, they developed a more efficient computer routine. Finally, Fish 
and Lee 4,26,27 and Botting, Vizzini and Lee 19 used a three-dimensional assumed stress 
hybrid finite element approach in order to obtain more accurate stress predictions than 
conventional displacement based elements, which are used in the previously mentioned 
studies. Extensive mesh refinement at the end of the ply drop-off was incorporated into 
their analyses. However, each element spanned the half-width of the quarter symmetry 
model and therefore the benefits of a three-dimensional model were not fully utilized. 

Plane Stress and Strain Finite Element Analyses 

Because of the computationally intensive nature of three-dimensional finite element 
models, most of the authors reduced the domain of the problem to two dimensions. 
Wisnom 17 was the only author who used plane stress finite elements in his study of 
unidirectional glass/epoxy laminates, which do not suffer from free edge effects. This 
approach predicts the stresses along the edge of a specimen, where a state of plane stress 
exists. Salpekar, Raju and O’Brien 28 employed the plane strain assumption in their two- 
dimensional analysis. 

Quasi Three-Dimensional Finite Element Analyses 

Kemp and Johnson , 14 Curry et al ., 20,21 Wu and Webber , 8 Wu , 9 Kubr , 7 and Messick 
et al . 2 all reduced the domain to a cross-sectional plane parallel to the load axis for 
transverse ply drops through the elimination of the dependence of the strains on one 
coordinate. Kemp and Johnson and also Kubr assumed a displacement field independent of 
the coordinate normal to the plane of analysis. Curry et al., Wu, and Wu Webber assumed 
the strains to be independent of the coordinate normal to the plane of the analysis, which 
restricts but does not preclude dependence of the displacements on that coordinate. Finally, 
Messick et al. used the assumption of axisymmetry in their analysis of the shuttle solid 
rocket booster casing. In each of these models, although the domain was reduced to two 
dimensions, the displacement normal to the plane of the model was still included and 
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therefore these models have five or six nonzero components of strain. This type of finite 
element analysis will be referred to as quasi three-dimensional in this literature review. 

There are several features of these quasi three-dimensional analyses that are worth 
mentioning. In order to arrive at a solution that took into account the out of plane boundary 
conditions imposed by an actual test fixture with simply supported edges, Curry et al. first 
used a plate finite element routine with a two step change in thickness to obtain boundary 
conditions for their quasi three-dimensional model of the ply drop-off region. In their 
analysis of external ply drop-offs, Wu and Webber included a thin resin layer between the 
continuous and dropped layers and then in his follow-up paper, Wu added nonlinear 
material response for the resin layer. Messick et al. used submodeling in their analysis of 
the composite shuttle solid rocket casing because of the large number of plies involved. 
Salpekar et al. as well as Wisnom refined their mesh extensively at the tips of the dropped 
plies in order to allow for accurate strain energy release rate calculations in their fracture 
analysis. Fish and Lee 4 used an assumed stress hybrid element approach to examine the 
free edge effects in a dropped ply specimen. 

The remainder of the studies incorporating quasi three-dimensional finite element 
analyses examined longitudinal ply drops. Because of the orientation of the ply drop in 
these cases, the plane of the finite element mesh was oriented transverse to the principle 
loading direction. Chan and Ochoa 1011 used displacement assumptions similar to Curry et 
al., allowing for extension in the direction of, and torsion about, the out-of-plane normal 
(the specimen longitudinal axis), as well as bending about the transverse axis (the specimen 
width-wise axis). Vizzini 13 performed a similar analysis in tension using assumed stress 
hybrid elements. 

1.1.3 Failure Analysis 

In order to predict the performance of the various laminates studied, some kind of 
failure analysis was applied and often more than one was used. Frequently, different 
criteria were used for prediction of in-plane and out-of-plane failure of the plies as well as 
for out-of-plane failure between plies (delamination). Most studies used some type of 
strength-of-materials approach which compares the local stress or strain state to material 
strength allowables. Others used criteria based on fracture mechanics where the strain 
energy release rate is compared to some toughness parameter to determine whether a flaw 
will grow and lead to failure. 
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Strength-of-Materials Approach 

Starting with the strength of materials approaches, several authors used some form of 
the maximum stress criterion. Daoust and Hoa 3 used it to predict the location of both in- 
plane and out-of-plane failure although no experiments were performed to verify the 
results. Fish and Lee 26 27 applied it in instances where one stress component was 
dominant. For a case where the interlaminar shear stress was dominant in the interply resin 
surrounding the ply drop-off, the criterion was applied in combination with interlaminar 
stress averaging. In comparison to experiments, Fish and Lee found that with a stress 
averaging distance of one ply thickness along the interface, consistently accurate 
delamination initiation predictions were obtained. Curry et al. 20,21 used the maximum 
stress criterion in the resin rich region ahead of the ply drop, but failure initiation was both 
predicted and experimentally found to occur at the interface between the dropped and 
continuous plies. In their purely analytical paper, Kemp and Johnson 14 used the maximum 
stress criterion in terms of the principal stresses in the resin rich region and applied it to the 
interply resin surrounding the dropped plies as well. Lagace and Cannon 15 - 16 applied the 
maximum stress criterion in their analysis; however, only in-plane stresses were 
considered, and as mentioned previously, the laminates they chose tended to be unaffected 
by the ply drop-off. In their analysis of the shuttle solid rocket booster casings, Messick et 
al. 2 found a good correlation between the magnitude of the peak value of the axial stress, 
which did not change location under minor design changes, and the failure load of the 
specimen. 

While they used the maximum stress failure criterion for resin failure, Kemp and 
Johnson used the Tsai-Wu criterion for intralaminar failure prediction. Curry et al. also 
used Tsai-Wu for intralaminar failure prediction but with a refinement developed by Hashin 
that eliminates the influence of a tensile strength parameter on a compressive failure mode 
and vice versa. In addition, Curry et al. used an interlaminar criterion based on matrix 
failure modes, also developed by Hashin. In combination with their quasi three- 
dimensional finite element analysis, these criteria underestimated experimentally determined 
failure loads by more than 30%. In addition to applying maximum shear stress to the 
interply resin. Fish and Lee also used a form of the Tsai-Wu criterion that considered only 
interlaminar stresses for out-of-plane failure of the plies as well as the interply resin. This 
criterion did not give good correlation with experiments for out-of-plane ply failure, but 
gave good results when applied to the interply resin with a longitudinal stress averaging 
distance of one ply thickness. Vizzini 13 used a combination of Tsai-Wu and von-Mises 
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criteria in his analysis of longitudinal discontinuities, including ply drops. The Tsai-Wu 
criterion was applied to the case of alteration of the ply orientation near the free-edge in a 
model without a resin zone at the discontinuity. The von-Mises criterion, which is an 
isotropic criterion, was applied within a resin zone for a modified model of the edge 
alteration as well as to the resin zone of the ply drop-off, and provided better strength 
prediction when compared to experiments. 

Wu 9 was the only author to use the maximum strain failure criterion in his study of 
external ply drop-offs. Coupled with his finite element analysis employing nonlinear 
material behavior, he obtained reasonable delamination prediction compared to 
experiments. 

Fracture Mechanics Approach 

Salpekar et al. 28 and Trethewey et al. 22 applied linear elastic fracture mechanics to 
predict failure. In both cases, their analyses were used to calculate the strain energy release 
rates for modes I and II which were compared to the corresponding fracture toughness of 
the matrix. However, neither of these papers included an experimental investigation and 
therefore the accuracy of these methods has not been verified. 

Wisnom 17 ’ 18 concludes the strain energy release rate is the critical parameter affecting 
delamination in dropped-ply unidirectional glass-epoxy specimens. In his experiments the 
onset of delamination was more difficult to measure objectively than the delamination 
propagation. The finite element analysis required considerable effort to compute initial 
strain energy release rate, and it was found that the initial strain energy release rate for the 
specimens with cut central plies was higher than the propagation value, while the opposite 
was true for the dropped ply specimens. (The propagation value corresponds to the 
uniform value of the strain energy release rate attained with increasing delamination crack 
length.) Moreover, a simple formula was derived to compute the applied stress for the 
delamination propagation which does not require a finite element analysis, but was 
confirmed by finite element analysis to be accurate. The experimental results for dropped- 
ply specimens correlated reasonably well using the simple formula to predict the applied 
stress for delamination propagation when the critical value of the strain energy release rate 
was measured from the straight specimens with cut central plies. For tension loading, the 
simple formula gave conservative predictions that were within 27% of measured values in 
Ref. 17, and within -20% to 3.9% in Ref. 18. For compression 18 , the simple formula 
gave conservative predictions that were within -5.5% to 1.9%. However, it was observed 
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that the critical strain energy release rate was not a material parameter since it increased as 
the thickness of the straight specimen with cut central plies increased. Also, tests on scaled 
specimens, in which the ratio of the number of cut central plies to continuous plies is 
constant, showed a scale effect since the delamination strength decreased with increasing 
thickness. (This phenomenon was also noted by Brewer and Lagace 29 in their study of 
delamination at straight free-edges.) Wisnom 18 investigated the use of the point stress 
criterion and the average stress criterion applied to the straight specimens with cut central 
plies, and found that the distance parameters in these criteria are not material constants 
either, but vary significantly with the specimen geometry. 

1.1.4 Summary Comments 

For studies in which results from analyses are compared to experimental results, 
deductions can be made regarding the merits of various approaches to analysis and failure 
prediction. A few features were common to two of the studies that obtained a good 
correlation to experiments. Fish and Lee 26 - 27 and Vizzini 13 used hybrid finite element 
analyses and applied their strength-of-materials failure criteria within interply resin zones. 
Fish and Lee found that better prediction of delamination was made when they applied their 
Tsai-Wu based criterion to an interply resin zone that was 10% the thickness of a lamina 
than when applied within a lamina. They also incorporated a stress averaging distance of 
one lamina thickness. Vizzini applied the von-Mises criterion to the resin rich zone 
bordered by the ends of the dropped plies and the two continuous layers and did not use 
stress averaging. However, his analysis produced stresses at only one point for each 
element which probably had the effect of averaging. Both of these analyses also employed 
finite element meshes that contain step changes in slope of the lamina which is not 
consistent with the actual geometry of a ply drop-off. In a follow up of his master’s 
thesis 20 , which had abrupt slope changes, Curry 21 found that smoothing the changes in 
slope did have a fairly large effect on the stress state. The interlaminar shear stress was 
reduced by about 30% through this smoothing in some cases. 

Wisnom 17 also was able to closely match his analysis with experimental results. His 
use of a strain energy release rate approach proved very accurate for laminates that were the 
same thickness as the one he used to measure the fracture toughness. However, he found a 
dependence of that toughness parameter on the thickness of the laminate. 

Other studies that attempted to predict failure with a comparison to experiments were 
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those by Curry, 20 - 21 Wu, 9 and Cannon. 15 - 16 Curry et al. underestimated the experimental 
failure load by 30%. Including a stress averaging distance or a characteristic distance would 
probably have improved on this result. Wu’s analytical prediction was off by about 25% of 
the failure load which he attributed to not including a fillet in the comer of the external ply 
drop-off. The analysis by Cannon did not attempt to model the geometry in the ply drop-off 
region accurately, and therefore the poor correlation with experiment in predicting failure is 
understandable. 

This review of the literature indicates that delamination is the predominant mode when 
dropped-ply laminates fail prematurely. Much has been written about delamination in 
laminated composite materials, mostly dealing with the firee-edge problem. In their paper 
presenting a quadratic stress criterion for predicting delamination initiation. Brewer and 
Lagace 29 described delamination as an initiation and growth process. After initiation, the 
delamination can undergo stable or unstable growth eventually leading to in-plane and final 
failure. While delamination initiation may not be the ultimate failure of a laminate, it is 
nevertheless important to predict. 

While the strain energy release rate approach has shown the ability to predict initiation 
and growth characteristics of delamination, there are major drawbacks to the method. In 
predicting initiation, for instance, the location of the delamination must be assumed in order 
to determine the strain energy release rate for a crack located there. Also, as mentioned 
previously, the critical strain energy release rate appears to be a function of laminate 
thickness. Finally, O’Brien 30 found that the critical strain energy release rate is a function 
of the percentage of Mode I present which requires a detailed analysis such as finite 
elements to determine. 

The main advantage of the strength of materials approaches is their ease of application, 
thus allowing the application of the failure criteria to more of the laminate at minimal 
expense. This advantage makes locating the point of delamination initiation more efficient 
than with the strain energy release rate approach. However, as stated previously, 
delaminations can undergo stable growth and mechanics of materials approaches cannot 
discern this characteristic of a delamination. Therefore, both of these approaches have then- 
uses in failure prediction. 

1.2 Objective 

The objective of this research is to study the use of an assumed stress, or stress-based, 
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method of approximation to predict interlaminar stress distributions in the vicinity of the ply 
termination. The aim is to develop a model accurately reflecting the geometry and material 
discontinuity of the dropped ply laminate such that distributions of the interlaminar stresses 
may be used with confidence in a point stress criterion or an average stress criterion. Stress 
based methods are discussed in Chapter 2, and Chapters 3 and 4 present specific models 
based on the method of analysis chosen in Chapter 2. The accuracy of the numerical 
methods used to solve the models is discussed in Chapter 5. Parametric studies are given in 
Chapter 6, and concluding remarks are given in Chapter 7. 
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Chapter 2: Choice of Method of 
Approximation 


This chapter begins by presenting a discussion of various approaches that were 
considered and the reasons for the choice of analysis that was made. Then, a general 
outline of the approach is presented. Following that, specific formulations are detailed that 
are of interest for the analysis of dropped-ply laminates. 

2.1 Methods of Approximation 

The most accurate way of solving for the static linear response of a solid body is to 
develop an exact analytical solution to the governing partial differential equations of three- 
dimensional linear elasticity satisfying the specified boundary conditions. However, for 
most problems, including laminates containing ply drop-offs, exact solutions are intractable 
and some form of approximation must be used. There are two general types of methods for 
approximating the continuous response of structures that are too complicated for an exact 
solution. The first involves developing an approximate solution to the governing 
differential equations of elasticity. The second uses the variational principles of solid 
mechanics and assumptions regarding the stress and or displacement fields to either derive 
the corresponding governing Euler equations or to arrive at the approximate solution 
directly. 

2.1.1 Approximate Solution to the Governing PDEs of Elasticity 

The primary approach to solving the governing equations of elasticity that is suitable for 
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complicated geometries such as ply drop-offs is the finite difference method. The finite 
difference method approximates the field variables on a mesh of points distributed 
throughout the domain. Approximations of the partial derivatives in terms of the field 
variables at the grid points, known as difference quotients, are substituted into the PDEs 
giving the difference form of the PDEs. These difference equations, applied over the 
domain, along with discretized forms of the boundary conditions, give a system of 
algebraic equations that can be solved for the grid point values of the field variables. 

There are several drawbacks to the finite difference method including inaccuracy of the 
derivatives of the approximated solution, difficulty in imposing boundary conditions on 
curved boundaries, and the inability to employ nonuniform and nonrectangular meshes. 31 

2.1.2 Variational Methods 

Variational methods are much more common in dealing with complex structures than 
methods that operate direcdy on the governing equations of elasticity. Based on the 
variational principles of solid mechanics, variational methods can be applied in many ways 
and can be tailored to the problem at hand and the results desired. 

One way in which the variational principles may be applied is to develop a structural 
theory based on some type of assumption regarding the field variables. An example of this 
approach is the development of classical beam theory through application of the principle of 
minimum potential energy. In this instance, the displacement field is restricted to a form 
explicit in the cross-sectional coordinates and substituted into the principle, which is then 
integrated with respect to the cross-sectional coordinates, giving the governing equation of 
classical beam theory as the Euler equation. A solution to a problem can then be obtained 
through solution of the boundary value problem consisting of the governing equations and 
the associated boundary conditions of the structural theory. 

Another class of variational methods, which obtain approximate solutions directly from 
the variational statement, are known as direct methods. They include the Ritz method and 
also the weighted residual methods of Galerkin, Kantorovich and Trefftz. It should be 
noted that the weighted residual methods can also be applied in instances where a 
variational statement does not exist and are therefore not exclusively direct methods. 
Because the Ritz method is applicable to any problem for which a variational statement 
exists and gives results that are in general superior to the weighted residual methods 32 , it is 
far more common in structural mechanics. The most popular use of Ritz type methods is in 
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the development of finite element approximations. 

In the finite element approach, the domain is divided into subdomains or elements with 
the form of the solution assumed within each element in terms of discrete values and 
polynomial interpolation functions. Application of the Ritz method to the appropriate 
variational principle in a piecewise manner results in a linear system of equations that is 
then solved for the discrete unknowns. While they may at first seem similar, these methods 
are generally recognized as being superior to finite differences. 

2.1.2.1 Conventional Finite Element Methods 

Within the field of finite elements as applied to solid mechanics there are several types 
of models which differ in the form of variational principle employed and the field variables 
that are assumed . 33 34 Using the terminology of Pian and Tong , 33 those based on the 
conventional variational principles of minimum potential energy, minimum complementary 
energy and the Hellinger-Reissner principle are the compatible, equilibrium and mixed 
models respectively. 

In the compatible models, the assumed field variables are the displacements which are 
required to be continuous within and between elements and therefore compatible (or 
kinematically admissible). They are the most commonly used of the finite element models 
because of their inherent ease of development for most applications and efficiency of 
computation. Their ease of development is due to the relatively loose restriction of 
continuity on the assumed displacement field. However, this simplicity does not come 
without a cost, namely a loss of accuracy in predicting stresses. This loss of accuracy is 
due to the fact that equilibrium of the stresses within the elements is satisfied only in an 
integral sense as is traction continuity between elements. 

The equilibrium models are derived from an assumed equilibrating stress field that must 
also satisfy interelement traction equilibrium, i.e., a statically admissible stress field. These 
models tend to have improved accuracy for stress prediction as compared to compatible 
models. However, satisfying the differential equations of equilibrium within an element 
and traction continuity between adjacent elements is very difficult for elements of varying 
shapes. Another drawback is that displacement continuity between adjacent elements is 
satisfied only in an integral sense. 

The mixed models assume both the stress and displacement fields with no a priori 
equilibrium requirement for the stresses nor compatibility requirements for the 
displacements. Conditions on the stresses and displacements between adjacent elements. 
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which are developed from the condition that the functional must be uniquely integrable, are 
also more flexible. Eight possible combinations result from the requirement that either the 
displacement or traction must be continuous in the normal and two tangential directions 
relative to the interface . 33 This feature makes the mixed finite element models useful for the 
analysis of plates and shells where difficulties in developing interpolation functions that 
meet continuity requirements on the derivatives of the out-of-plane displacement can be 
avoided. Moreover, it also provides benefits in stress analysis by allowing assumption of 
stresses, which according to Pian and Tong 33 can improve their accuracy, without the 
complication of having to satisfy equilibrium a priori. 

2.1.2.2 Hybrid Finite Element Methods 

Beyond the conventional variational principles are the modified variational principles 
which allow for relaxed continuity requirements between elements. The models based on 
these principles are the hybrid models and are distinguished by the independent assumption 
of field variables within the elements and along the element boundaries which are allowed 
by the modified principles . 35 - 36 Derived from the conventional variational principles 
through the use of Lagrange multiplier terms that apply the conditions on the tractions and 
displacements between elements as constraints, the modified principles of potential energy 
and complementary energy lead to the hybrid-displacement and hybrid-stress models 33 
while the modified Hellinger-Reissner principle leads to the hybrid-mixed models . 34 - 37 In 
the hybrid-displacement models, the displacement field is assumed within the elements 
while the tractions and also possibly a separate form of the displacements are assumed on 
the element boundaries. The hybrid-stress model assumes interior stresses and boundary 
displacements. The hybrid-mixed models use assumed element stresses and displacements 
with separate boundary displacements and possibly tractions. 

The hybrid-displacement models take three basic forms, two based on the assumption 
of element displacements and boundary tractions and one that assumes element 
displacements and both boundary tractions and displacements . 34 - 38 The difference between 
the first two forms is that one bases its displacement assumptions on nodal displacements 
which guarantees compatibility with neighboring elements while the other uses general 
displacement parameters. The third form is a combination of the first two in that it assumes 
the element displacements in terms of general displacement parameters while the boundary 
displacements are in terms of nodal displacements. These models offer flexibility in dealing 
with interelement displacements and like the conventional mixed models, are typically 
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applied to plate and shell problems. 

First introduced in 1964 by Pian 39 , hybrid-stress models are typically applied where 
more accurate stresses are desired. As in the equilibrium model, this hybrid model uses an 
assumed equilibrating stress field within the elements which enhances stress accuracy. 
However, it avoids the problem of meeting interelement traction equilibrium through the 
modification of the variational principle. Instead, the boundary displacements are assumed 
in terms of nodal values such that they satisfy interelement continuity. Of the hybrid 
models, these have been applied to the broadest range of problems, from multilayer plates 
and shells to fracture mechanics. 

The hybrid-mixed models, like the hybrid-displacement models, have been applied to a 
relatively small number of problems as compared to the hybrid-stress models. In fact, it is 
considered to be impractical in the field of linear solid mechanics to assume all of the 
variables allowed by this approach, i.e. separate internal and boundary values of the 
stresses and displacements. 34 - 36 Instead, this approach is recommended for dynamic 
problems and incremental analysis of finite deformation problems. 

Application of these finite element models to the analysis of laminates can be achieved 
in two ways. The first is to simply apply elements developed for the analysis of continuous 
media, i.e. continuum elements, to the laminate. In this approach, the individual plies can 
be modeled by one or more elements through the thickness, or groups of plies can be 
modeled by one element through their collective thickness using some type of smearing 
scheme to arrive at an approximation of their constitutive behavior. The second approach is 
to develop a laminate finite element that takes into account the multilayer nature when 
developing the element approximation of the field variables. 

A few studies have been done regarding hybrid-stress or mixed finite element analysis 
of laminated plates that use the second approach rather than simply discretizing the layers 
with ordinary continuum elements. Spilker 4 * 3,4 * formulated hybrid-stress laminate elements 
that are isoparametric in the transverse coordinates with stress fields that equilibrate within 
each layer and satisfy traction continuity between adjacent layers. However, this approach 
does not allow for tapered layers as is required in analyzing dropped-ply laminates. Cook 42 
also developed hybrid-stress elements, but they were based on plate theory and also cannot 
account for taper. Most recently, Shi and Chen 43 developed a global-local model that uses 
compatible elements in the global portion and hybrid-stress elements in the local region. 
Once again, however, taper cannot be accounted for in their model. 

The difficulty in formulating a hybrid-stress model of a tapered laminate using this type 
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of approach lies in satisfying both the equilibrium equations on the interior of each layer 
and traction continuity at interfaces for layers with curved boundaries. For a finite element 
approach, one of the two must be relaxed, which is essentially the same as going back to 
treating the laminate like any other domain and using standard continuum elements. 

2.1.3 Criteria for Interlaminar Stress Approximation 

For accurate prediction of interlaminar stresses from an approximate continuum model, 
the method of approximation should meet certain requirements. Specifically, the stress field 
should satisfy equilibrium within the layers of the laminate as well as equilibrium of 
tractions, also referred to as traction continuity, between layers. In addition, it would be 
desirable to satisfy displacement continuity between layers. 

Of the finite element methods presented, those that come closest to meeting these goals 
are the stress based methods of equilibrium and hybrid-stress and also the mixed methods. 
While the equilibrium models meet the requirements regarding the stress assumptions, 
development of such elements is difficult and displacement continuity between layers is not 
met. The hybrid-stress models do not enforce traction continuity and developing a mixed 
model that meets these requirements runs into the same difficulties as the equilibrium 
models. Therefore, it is felt that a different approach is necessary, which is the subject of 
this dissertation. 

2.2 Approach 

The approach chosen is based on work by Pagano 44 ’ 45 and develops a structural theory 
in a similar manner to the development of classical beam theory using the principle of 
minimum potential energy. However, instead of assumptions on the displacements, the 
present theory makes assumptions on the stresses. In addition, the Hellinger-Reissner 
stationary variational principle will be used instead of the principle of minimum potential 
energy 

Development of the theory begins with division of the domain into layers with the 
stress field assumed within each layer. The dependence of the stresses on the thickness 
coordinate is expressed explicitly with considerations for the satisfaction of equilibrium and 
traction continuity between layers. The Euler equations are obtained from the stationarity 
condition with traction and displacement continuity imposed. Solution of the governing 
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equations is by a one step finite difference scheme. Pagano’s presentation of this type of 
approach was first developed for analyzing flat laminates. 44 Later, he modified the 
approach to handle layers defined by curved meridians in axisymmetric bodies 45 


2.2.1 Variational Principle 

The starting point for the derivation of this analysis is the Hellinger-Reissner functional 
for a domain composed of N subvolumes and no body forces. In cartesian coordinates. 


n, = £j v [<W«d-«H)P dV * -I w«-J s -W* < 2 -'> 
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in which the displacements are denoted m,, the strains £ i; , the stresses cr, the surface 
tractions T jt and the complementary strain energy density W. Repeated indices are summed 
from one to three in the usual indicial notation and the tilde symbol (~) indicates a 
prescribed value. The portion of the surface with tractions prescribed is denoted by S a 
while S u represents the portion with displacements prescribed. In this form, the functional 
assumes that the displacements are continuous between subvolumes, i.e., the strains in the 
first term of the volume integral can be represented without the use of the delta function. 
Without this condition, an additional integral along the border between subvolumes must be 
included in order to account for the presence of the delta function in the strain field.-^ In 
addition, the tractions would have to be continuous for the functional to exist if the 
displacements are not continuous. 

Stationarity of the functional with respect to all possible variations in the displacements 
and stresses gives the governing equations of elasticity as its Euler equations. The 
vanishing of the first variation is given by 
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for all admissible 5o tJ and 8u r Beyond symmetry of the stress tensor which is required in 
order to arrive at the governing equations of elasticity from the principle, admissibility of 
the stresses and displacements is governed by the boundary conditions and the unique 
definability of the functional of Eq. (2-1). As mentioned previously, for the functional to be 
uniquely defined, the 0 {j £ tJ term must be integrable throughout the domain. Continuity of 
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displacements throughout the domain ensures this without the need for additional terms in 
the functional. Therefore, admissible variations in the displacements are to be continuous 
within subvolumes and at interfaces and vanish on 5„. As stated in the objective, it is 
desired to develop a theory in which tractions are continuous between layers. Therefore, 
admissible stress variations are to be symmetric (5o,j = &7 /( ), continuous in subvolumes, 
satisfy traction continuity at interfaces and vanish on 5^ 

In the cartesian reference, the variation of the strains written in terms of the 
displacements is 

Se u =i( Su i.j + Su j.i) ( 2 - 3 ) 

where the comma convention is used to indicate partial differentiation. Therefore, by 
substituting (2-3) and incorporating symmetry of the stress tensor: 

l °„K dV , = 1, A. 0-“) 

Applying a rearranged version of the product rule of differentiation and then using Gauss’ 
divergence theorem and Cauchy’s relation for surface tractions in terms of stresses leads to 
the result: 



where S k represents the surface of each subvolume V* and i; = ojij. 
Substituting (2-5) back into (2-2) yields 
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The new summation term results from the part of the surface integral in Eq. (2-5) taken 
over the portion of S k , which is denoted by /,, that borders another subvolume where 
neither tractions nor displacements are prescribed. In this term the value of M will depend 
on the geometry of the domain chosen and the + and - superscripts denote the two sides of 
the interface. 

From the previously stated condition that the displacement field is continuous between 
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subvolumes, the interface displacement variations are no longer independent, i.e. 

5u+ = 8 u~, and the internal interface term becomes the condition of traction continuity. For 
a domain composed of layers, there will be one less interface than there are layers. 
Numbering these layers sequentially such that layer k is adjacent to layer k + 1 for all k and 
denoting the interface between these adjacent layers by I k , the principle takes the final form 
for application to laminates. Thus, 
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2.2.2 Assumed Stresses 

Development of the assumed stress field is outlined for the general case of an 
orthogonal coordinate system with dependence of the stresses on all three coordinates. It 
should be noted here that while this procedure is applicable to all such coordinate systems, 
it becomes very involved for most cases. Use of this more general form is simply for 
illustrative purposes and the models developed subsequently are for specific applications of 
the procedure in simple coordinate systems. 

For the general orthogonal coordinate system (a p a 2 ,£), it is required that the surfaces 
that define the boundaries of the layers be expressed as £ = g{a v a 2 ) where g is a smooth 
function. The lateral boundaries would be defined by constant values of either a, or a 2 
(see Fig. 2.1). A multilayered domain could therefore be defined by a family of surfaces in 
the thickness direction and by constant values of ct l or cc 2 at the lateral edges. 

Considering a single layer bounded by £ = C l (cc l ,cc 2 ) and £ = C 2 (® v a 2 ^ £2 > Ci 

for all ctj and cc 2 within the domain, the stress field within that layer is assumed in the form 

a > (a 1 ,a 2 ,C) = CTj(a 1 ,a 2 )/j i) (a„a 2 ,C). nosumoni (2-8) 

where o i are the six stress components in contracted notation, o u are functions only of a, 
and a 2 to be determined through application of the principle, and fj‘ ] are shape functions 
with explicit dependence on £ and known implicit dependence on a, and a 2 through 
functions £, and £ 2 . The requirements that 
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<j i (a ] ,a 2 ,C i ) = cr i] (a l ,a 2 ) 
cr,(a 1 ,a 2 ,C 2 ) = <r i2 (a 1 ,a 2 ) 

arc imposed on the stress assumptions in order to facilitate interlayer traction continuity. 



Fig. 2.1 Section of a uniform thickness layer in (a lt a 2 ,C) orthogonal coordinates. 

In arriving at the specific form of the shape functions, the “in-plane” stresses <J U , <r 22 
and <t 12 are assumed to have a specific dependence on £. Based on the elasticity equilibrium 
equations, the through the thickness distribution of the remaining stress components are 
determined. This procedure involves substitution of cr n , <r 22 and <r 12 and integrating with 
respect to £ in order to determine the dependence on £ that the remaining stresses must have 
if non-trivial values of the stress variables o u are to satisfy equilibrium. 

2.2.3 Governing Equations 

The form of the complementary strain energy density W for linear-elastic materials in 
terms of the stresses in contracted notation is 

W = {S IJ G I C J (2-10) 

where S l} are the compliance coefficients. This expression is substituted into Eq. (2-7) 
along with the assumed stress field (2-8) and its variation. 
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( 2 - 11 ) 


8o ( ft , ft 2 , — S<Jj (ftj , cc 2 )fj (®i ’ C) 

and Cauchy’s relation for the surface tractions in terms of the stresses. Integration with 
respect to £ is then carried out. Because the displacements have unknown dependence on £, 
integrals through a layer of the displacements weighted by functions of £ are defined as 
new unknowns in the model. 

The Euler equations of the resulting form of the principle are one of two types. Those 
associated with the variation in < j u are the compatibility/constitutive equations, which may 
contain derivatives of the weighted displacements. Those associated with the variations in 
the weighted displacements are the equilibrium equations, which may contain derivatives of 
the stress variables. It should be noted that, while satisfaction of the equilibrium Euler 
equations implies satisfaction of the elasticity equilibrium equations within the layers, 
satisfaction of the Euler compatibility/constitutive equations only implies satisfaction of the 
corresponding elasticity equations in an integral sense within the layers. The system of 
Euler equations contains both first order linear ordinary differential equations and linear 
algebraic equations. A system containing both differential and algebraic equations is 
referred to in the literature as a system of DAEs. The solution of the systems will be 
detailed following the development of specific models using this approach. 
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Chapter 3: Generalized Plane 
Deformation Model 


The first version of the model is developed for problems of generalized plane 
deformation in cartesian coordinates. Therefore, (a p a 2 ,^) = (x,y,z) for this formulation. 
In this class of problems, the loading, geometric and material properties, and hence stresses 
and strains are independent of y. The solution domain is thus reduced to the x-z plane 
where the x-axis is the longitudinal direction and the z-axis is the thickness direction. These 
problems differ from generalized plane strain in that bending about the x- and z-axes and 
torsion about the y-axis are allowed. A generic dropped ply laminate represented in a 
generalized plane deformation model is shown in Fig. 3.1. In this figure P denotes the 
extensional force along the y-axis, M x is the bending moment along the x-axis, M z is the 
bending moment along the z-axis, and T denotes the torque along the y-axis. Actions P, 

M x , M z , and T are independent of the y-coordinate, and they lead to displacements that 
depend on coordinate y in an explicit manner. 

3.1 Displacement Field 

For this class of problems, the most general form of the displacement field (see page 
104 of Ref. 46) is 

u(x, y, z) = - 1/2 Ay 2 + Oyz + U( x, z) + u' 

v(x, y, z) = y ( Ax + Bz + C) + V{x, z) + v' (3- 1 ) 

w{x,y,z ) = -1/2 By 1 - 6xy+ W/x,z) + vi/ 
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where A and B are the negative bending curvatures in the x-y and y-z planes respectively, 
C is the normal strain at x = 0 and z = 0 , and the product 9y is the rotation of a cross 
section about the y-axis. A, B, C, and 6 are all constant over the domain and are the 
prescribed out-of-plane deformation quantities. Quantities A, B, C, and 0 can be related 
through Hooke’s law and static equivalence to actions P, M x , M z , and T of Fig. 3. 1. The 
functions U, V and W are the unknown portions of the displacements which are functions 
only of x and z, and u', v', and w' are the rigid body displacements composed of three 
translation and three rotation modes. 



Fig. 3.1 Generalized plane deformation. 


It is worth noting that this displacement field does not include all possible out-of-plane 
deformations (those whose displacements involve dependence on y) which have strains that 
are independent of y. For instance, a uniform state of simple shear that would correspond 
to a term consisting of the product of a constant times y in either the u or w displacement is 
not included. Because the tractions on the edges parallel to they axis are arbitrary at this 
point, such a deformation is precluded by the fact that it would produce a stress field that 
varies with y. 

Based on the displacement field in Eqs. (3-1), the engineering strains, which appear in 
the first term in the volume integral of Eq. (2-7) when the summation over i and j is earned 
out, have the following forms: 
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£ x = u x = U x 

£ y = v y = Ax + Bz + C 
£ z = ", = W t 

r yl = v z + w y = V z -6x 

y xz = u .z + ", = u . z + w x 
Y» = “y + v x =V x + dz 

Note that these representations of the strains are independent of y as required. 


3.2 Assumed Stress Field 


A generic layer for the generalized plane deformation model is shown in Fig. 3.2. Eq. 
(2-8) for this problem has the form 

o i (x,z) = (Jjix )fj‘\x,z) , no sum on i (3-3) 


As outlined in Sec. 2.2.2, the in-plane stresses <7„, cr^ and (Gi, <J 2 and cr 6 in contracted 
notation) are assumed to have linear dependence on z. Based on the y-independent form of 
the elasticity equations of equilibrium without body forces. 


°xx,x 


G 

+ d„ =0 

xy,x 

yz,z 

G 

X Z,x 

+ <r«, = ° 


(3-4) 


the through the thickness distribution of the remaining stress components are determined. 
The through-the-thickness shear stresses o yz and o xz (<r 4 and <r 5 ) are found to have 
quadratic dependence on z and the thickness normal stress a a (<x 3 ) has a cubic distribution. 
In their most general form, the stresses can be expressed as 


0 2 (x,z) = a l (x) + b l (x)z 

<j 2 (x,z) = a 2 (x) + b 2 (x)z 

<x 3 (x,z) = a } (x) + b 2 (x)z + c i (x)z 2 +d i (x)z i 

cr 4 (x,z) = a 4 (x) + b 4 (x)z + c 4 (x)z 2 

<j 5 ( x , z) = a 5 (x) + b 5 (x) z + c 5 (x)z 2 

<J 6 (x,z) = a 6 (x) + b 6 (x)z 
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To arrive at stress assumptions of the form given in Eq. (3-3), the conditions of Eqs. (2-9) 
are used in the form 


<7,0 c,z,) = <T a (Jf) 

a,(x,z 2 ) = ct i2 (x) 



Fig. 3.2 Schematic of a typical layer. 

For the in-plane stresses Oj, o 2 and o 6 , these conditions lead in a straightforward 
manner to the following expressions: 


<t 2 


= <7, 



+ <7, 


= (7- 


2 , “2 


( Z-Zi y 
\ Z 1 ~ Z 1 J 
( 


12 


21 


l+<7. 


2 “ 2 , 


. 2 9 2 


22 




\ Z 1~ Z \. 


^6 ~ °61 


( z 2 - z ^ 

\ Z 2 ~ Z lJ 


+ o, 


62 


Z~Zj 

-2, 


i J 


(3-7) 


The out-of-plane stresses, on the other hand, are not uniquely determined by the two 
conditions of Eq. (3-6). In fact, there are an infinite number of possible choices for them. 
A simple method for determining a form for these stresses is to eliminate two of the 
coefficient functions from the general form for the stresses of Eq. (3-5) using the 
conditions of Eq. (3-6). Using this approach, the resulting form of the assumed stresses 
depend on which of the coefficient functions a„ b, etc. are eliminated through imposing 
Eqs. (3-6). By eliminating a 4 and b 4 from the <7 4 expression of Eq. (3-5), the result for <7 4 
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IS 




Z 2 

K z l — -M 


l+<X 


Z-Zy 


42 


V>2 


“2, 


+ c 4 (z 2 -z(z, + z 2 ) + z,z 2 ) 


1 J 


(3-8) 


The other two possibilities, elimination of a 4 and c 4 or b 4 and c 4 , result in quadratic 
polynomials in z as the shape functions multiplying <j 41 and <J 42 . The form given in Eq. (3- 
8 ) will be used because it contains the same linear shape functions as the in-plane stresses 
of Eq. (3-7). At this point, c 4 could simply be replaced by <r 43 and Eq. (3-8) would have 
the desired form given in Eq. (3-3). However, that would leave cr 43 without the dimensions 
of stress. In addition, the extremum value of (z 2 - z(z ( + z 2 ) + z,z 2 ) within the layer is 
— (z 2 - zj 2 / 4 which, for thin layers (z, = z 2 ), approaches zero. Both of these situations are 
avoided by replacing c 4 by cr 43 /(z 2 - z x ) 2 which is also a function only of x. The same 
approach is used for (T 5 , leading to 



r \ 

f \ 


z 2 -z 


Z"Z X 

p 

II 

£ 


+ ct 42 



\Z*i ~ Z\ j 


^Z 2 “ Zl J 

f \ 



z 2 -z 


r Z-Zj 

- °51 


+ <^52 



K z 2 ~ Z l j 


V Z 2 ~Z\ J 


+ o 


+ (X 


(z 2 -z(zi+z 2 )+z 1 z 2 ) 

(Z2-Zl) 2 

(z 2 -z(z i +z 2 ) + z ] z 2 ) 
(z 2 -z x ) 2 


(3-9) 


The shape functions associated with <t 45 and <J 53 are referred to as bubble functions because 
they are zero at z = z l and z = z 2 and therefore do not contribute to the surface values of the 
stresses. This designation is borrowed from finite element methods where a bubble 
function is one that vanishes on the boundary of the element. 

Using the same procedure on 0*3 that led to Eq. (3-8) gives the following expression: 


^3 “ ^31 



+ <x 


32 


/ \ 
z-z x 

^2 “ ^1 > 


+ c 3 (z 2 -z(z,+z 2 ) + z,z 2 ) 


(3-10) 


+4 3 (z 3 - z(z , 2 + Z[Z 2 + z 2 ) + z x z 2 {z x + z 2 )) 


Once again, this expression has the form of Eq. (3-3) by replacing c 3 by cr 33 and 4 3 by a M . 
However, the cubic shape function in Eq. (3-10), which is also a bubble function, has its 
third root at z = -(Zi + z 2 ). In order to separate its contribution to the distribution of <r 3 
within a layer from the quadratic term as much as possible, the cubic bubble function is 
modified to locate its third root at (z, + z 2 ) / 2 . This change also makes the cubic bubble 
function orthogonal to the quadratic one within the layer, that is, the inner product of the 
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two over the interval < z 5= z 2 equals zero. The scaling factor for this modified cubic 
bubble function that keeps the extremum value within the layer from approaching zero for 
thin layers is (z 2 - z,) 3 . Therefore, c 3 and d 2 are given by the following: 


<■3 



d 2 — 


2^34 

(*2-*,)* 


(3-11) 


The resulting expression for <r 3 is 


^3 ” ®31 
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+ °32 
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+ a. 


33 


(z 2 -z(z l + z 2 ) + z 1 z 2 ) 
(*2-Zl) 2 


+c. 
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(2z 3 - 3z 2 (z, + z 2 ) + z(zf + 4 z,z 2 + z 2 ) - z t z 2 (z, + z 2 )) 
(z 2 -Zi) 3 


(3-12) 


Therefore, the shape functions of Eq. (3-3) are given by 


yd) _ y( 2) _ y(3) _ y( 4) _ y(5) _ y(6) _ _£ £ 

z 2 — z t 

/d) _ A2) _ f(3) _ f (4) _ f (5) _ A6) _ Z~Zj 
h - J2 ~ h - h - h -J2 “ 

z 2 z L 

f(3) _ A 4) _ f(5) _ Z 2 -z(zi + Z 2 ) + Z|Z 2 
/3 “ h - h - I 

U 2 -Z 1 ) 

(3) _ 2z 3 - 3z 2 (z t + z 2 ) + z(z 2 + 4 Zl z 2 + z 2 ) - z,z 2 (zi + z 2 ) 

J4 — " “ ~ ™ 


(3-13) 


(«2 — Zl) 


These shape functions are used to define the assumed stress field for each layer of the 
domain with Zi and z 2 representing the lower and upper surfaces of a particular layer. It 
should be noted that the dependence of fj l) on x expressed in Eq. (3-3) is implicit through 
the presence of Zi and z 2 . More specific notation regarding the description of these surfaces 
is introduced in the following section. 
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3.3 Application of the Principle 


The governing equations for this model are developed through the procedure outlined in 
Sec. 2.2.3. Taking advantage of the independence of the integrands on y in Eq. (2-7) for 
generalized plane deformation, integration through a unit depth in y is carried out. This 
simply reduces the volume integral to an area integral in the x-z plane and the surface 
integrals to contour integrals in the same plane. Integration with respect to z is then carried 
out, reducing the solution domain to one dimension. Because of the dependence of Zi and 
Z 2 on x, Leibnitz’s theorem must be used on terms involving derivatives in x. 

At this point, the displacement terms U, V and W are the only remaining unknowns in 
the area integrals that are functions of z. In order to complete the integration with respect to 
z, the following definitions for weighted integrals are made within each layer: 

{g(x\g(x),g(x)) = p[h l ,h 2 ,h i ]g(x,z)dz (3- 14a) 

where 


f h = 

h = 


z 2 -z 

(z 2 -*i i ) 2 
z-z x 
(z 2 -Zi ) 2 

Z 2 -z(zi + 

(*2-*l) 3 


(3- 14b) 


and g represents a displacement component. The weighting functions h u h 2 and h 2 are 
simply the first three shape functions of Eq. (3-13) divided by the layer thickness (z 2 - z,). 
This scaling is done in order to retain the dimensions of length for the weighted 
displacements and to keep the integral of the weighting functions from approaching zero for 
thin layers. 

The weighting functions for the displacements in Eqs. (3- 14b) are different from 
Pagano’s 44 ’ 45 , which are simple powers of the thickness coordinate. This change was 
made to overcome problems with matrix conditioning encountered when numerically 
solving the governing equations for thin layers or layers located at moderately large values 
of z relative to the layer thickness. The reason for this change is illustrated in Fig. 3.3 
which shows plots of h x , h 2 and h 2 within a layer compared to Pagano’s 1, z, and z 2 . As 


Generalized Plane Deformation Model 


31 






layer thicknesses become smaller relative to their distance from the origin, of which Fig. 

3.3 is not an extreme case, z 2 approaches a straight line within the layer. Therefore, a linear 
combination of 1 and z can closely approximate z 1 within the layer. The choice of 
weighting functions given in Eqs. (3- 14b) clearly avoids this. 

To this point, the derivation of the governing equations has centered on a single layer. 
However, the final model is composed of several layers and requires a slight modification 
of the notation as well as the definition of a few additional terms. A schematic of an 
assembly of N layers is shown in Fig. 3.4 along with some new notation. The z-coordinate 
of the &th interface at the top of layer k is designated z*(jc) and the angle between the 
horizontal and the tangent to interface k is y k . From the geometry, z' k (x) = tan y k (x) , where 
a prime denotes an ordinary derivative with respect to x. The bottom of layer 1 and the top 
of layer N are part C s of the edge curve C, and ends jc = jc, and x-x 2 (x 2 > x, ) are part 
C of C (C= C +C e ). 



Fig. 3.4 Schematic of layer assemblage. 

Layer variables are denoted by the layer number in parentheses as a superscript. 
Additional subscripts 1 and 2 are used to denote dependent variables evaluated at the 
bottom and top, respectively, of a layer; e.g., u a> (x, z*_,(x)) = uf k> (x) and 
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u ik) (x,z k (x )) = u 2 \x ) . The resulting form of the variational principle for the structural 
approximation is 


f gift- + Xuf'SaV ~{G,SU + G,SU + G,SV + G,SV + G,SW + G S 5W 

+ r„ - xj"5u"' + (r„ - y„ 

+ [(r„ - iJ"’Sur +(r,,- f*f - ^r^'jsec 

-([(«, -s,)"’ 5T , ,v+(v 1 -v 1 ) m 5r; 1 >+(»' 1 *5!’]»cr. 

+[(“! +(v 1 -v 2 r5z;»> + ( W! -w 2 r<’]s e or„) (3 


+X[(og’ - <*; - <" + <C' W + (<C - <*; - °r ’ + )K' 

*=1 

+« - - < +i> + < +l) z*)^r]W 


(3-15) 


4{([K - *nf ^ +K - d n fSU^ + (a 6l - o a ) w SV*' 

+fo 2 5V ( *> +K -d 5l fsw^ +(o 52 -a 52 fSW ik) 

+K - -z^)) Q -([(^-^r > ^n ) +(0-d) m So&' 

+{v - P) W Sag +(v- vf k) 8&£ + (W- W) w S < + (w - #) (t> 
+(w-w) < * > 5<](z*-z t _,)^| =0 

Before describing the contents of Eq. (3-15), some comments on the terminology is 
needed. In the form presented here, the principle does not provide all of the Euler equations 
because of interdependence among some of the varied quantities. In the context of the 
variational principle, the word “term” will refer to quantities that may require additional 
consideration in that regard. The word “equation” will be reserved for the Euler equations, 
that is those associated with independent variations. 

The integral on jc from to x 2 in Eq. (3-15) contains compatibility terms associated 
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with the sixteen variations 5o u , equilibrium equations associated with the seven variations 
in weighted displacements, traction or displacement prescribed terms on C\ and internal 
interface condition equations in which Cauchy’s formula was used for the surface tractions. 
The remaining terms are boundary terms evaluated on C e of C. 

The compatibility terms in the variational principle of Eq. (3-15) contain the quantities 
fa, which in turn contain displacements evaluated at the interfaces, and Xu* which relate the 
weighted displacement quantities to the stresses. The nonzero fa are 


Pu = z ' k . , u\ k) 

Pn =-z>2* > 

< = 


= “ v i < * ) 

= V 2 ] 

^=-u[ k) + z k _ x w[ k) 

Hs2 =u[ k) -z k w ( j 

T 

II 

5L 

TS 

s ^ 

II 

1 

A- 


The Xu are 

jtf’-C + l ?’-■©<’ (3-17) 

in which weighted integrals of the compliances are defined by 

= (3-18) 

and the out-of-plane prescribed strains are contained in the quantities 

r® = l/6(z, - z._, )[3( At + C) + B( 2z k _, + z t )] 
r'« = V6(z. - Zt-i )[3(Ax + C) + B( z,_, + 2z„ )] 

r®=-V2(z,-z,.,)ex (3-19) 

C=l/6(z,-J,.,)( 2 z.-, + ^)« 

n 2 ‘’ = l/6(z,-z 1 .,)(z,_ 1 + 2 Zl )9 

Finally, the nonzero strain measures in terms of the weighted displacement quantities are 
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The equilibrium equations are simply G- k) = 0 where 


G?> 

G (k) 

G (k) 

<%' 

Gt' 

Gt } 

Gj k) 


= {z k - 

“ Zk- 1, 

X' 

+Z*-i( 

'o (k) - 
,11 

On)~ 

„ 4- n ,< * ) 

°51 ^°52 

_ a <*> 

°53 


= (z k - 

“ Zk-\ , 

X' 

+z*'« > -c 



X 


1 

II 

" 1 

K )f 

U 41 

+ <J tk) 

T U 42 

-< 

+c iK’ 

-<) 


= {z k ■ 

“ £*-l 

nr 

U 41 

u 42 

+ <C 

+z;(<c- 

<) 


= (z* ■ 

“Vi 

nr 

U 31 



+<45 >+ c 


< + <) 

= (z k ■ 

~ £fc-l 

nr 

— /T^) 

a 31 


+X 

+ < ) + z*'( 


^’-0 


= (z, -z M )oS , '-K>oS’ -2fe - C.X 


(3-20) 


( 3 - 21 ) 


3.4 Interlayer Traction and Displacement 
Continuity 

As mentioned previously, before the Euler equations of the principle can be determined, 
the interdependence among the variations of the unknowns must be accounted for. 
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Therefore, the conditions of admissibility stated in Sec. 2.2.1 on the stresses and 
displacements and their variations across subvolume (layer) interfaces are now imposed. 
Interlayer continuity of displacements is stated as 


u< k) =ut i+1) 
v<*> = v[ Mi 
w ( 2 k) = w[ k+l) 


(3-22) 


with the same conditions on their variations, which have already been accounted for. 
Traction continuity at interface z k is 

x - direction: (og* - ) = (og +1) - z*'< +1> ) 

y ~ direction: - z k O ( 6 k) ) = (< +1) - z*'< +1) ) (3-23) 

z - direction: (og J - z k <r£) = (<x* +1) - z t 'c7^ +1) ) 


It is assumed that the variations of the stresses satisfy traction continuity at the 
interfaces. Substitution of the variations of Eqs. (3-23) into the principle combines some of 
the compatibility terms, reducing their number by 3(2V— 1). Applying the continuity of 
displacements cancels the flu terms for the internal interfaces, eliminating the internal 
interface displacements from the formulation. The resulting interface compatibility 
equations are 


Xn+Xn +i) = 0; / = 1,3, 4, 5, 6 
^?+^a , +W)V» ) -0 ► 

Xa )+ z' k X ( 42=0 


k = l^>N-l 


(3-24) 


As an example of the procedure used to arrive at the surface condition equations, the y- 
direction conditions on the upper surface are developed. The two possible cases are 
v™ = v ( 2 N) or t'j } = 7 * 2 For the first case, and So ( 42 ] are independent giving the 

conditions 


XiV-z>\ m = o (3-25) 

;C+<> = o 

from which v[ N) can be eliminated and with minor manipulation leaves 
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(3-26) 


^+ 4 ^ = 0 

X*l + ^2 N) ~ 0 

For the traction prescription case where ' = T { y N 2 ’ , then - z^Sg^' = 0 (cos y N * 0) 
and the surface conditions are 

# 62 ) + 4 * 42 ) = 0 

o% } cosy w -<T^ ) siny iV = ^ > 

The first equation of Eqs. (3-26) and (3-27) are the same and therefore hold for either case 
while the second changes depending on the prescribed surface quantity. A similar 
procedure is used to arrive at the conditions in the x and z directions, however they are 
complicated somewhat by coupling introduced by the appearance of 1 and its variation 
in both directions for z' N * 0. The resulting conditions for the upper surface are 



(3-28) 

X^+z^+S^ = 0 

(3-29a) 

or 


cos y N - 0-j‘f sin y N = 

(3-29b) 

X^ + = 0 

(3 -30a) 

or 


($? cos y N - sin y N = fjf 

(3-30b) 


(3-3 la) 

or 


cos y N - sin y N = 1 

(3-3 lb) 


The conditions on the lower surface are similar and are presented later. 

The remaining compatibility terms are unaffected by the interdependence among the 
stress variables between layers. These terms lead to the layer compatibility equations which 
are defined for all N layers: 
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= Xn = =x£= X% = x£ = 0 (3-32) 

This completes the development of the Euler equations. 

3.5 Differential-Algebraic Equations Systems 

At this point, it is necessary to discuss the nature of this system of equations. While 
some of the equations contain derivatives of the unknowns and are therefore differential, 
others do not and are algebraic. Among the designations for this type of system are 
singular, constrained or descriptor. It will be referred to it as a system of DAEs 
(differential-algebraic equations) in the present work as that is most widely used in the 
contemporary literature on their solution. 

Solution of a system of DAEs may be thought of as integration of the differential 
equations and determination of the integration constants through application of the 
boundary condition equations. For a system of first order linear ODEs in which all 
unknowns appear as a differential, the number of these integration constants introduced 
corresponds to the order of the system. However, other factors come into play when 
dealing with DAEs and definition of an order for such a system becomes considerably more 
complicated. For instance, referring back to the above mentioned system of ODEs, which 
may be written in matrix form as 

Ay' + By = f , (3-33) 

consider the case for which some unknowns are present only as algebraic variables. 
Logically, if there are more differential equations than differential unknowns, A will have 
columns of zeros associated with the algebraic unknowns and will therefore be singular. 
Through elementary row operations on A, some of the rows can be driven to all zeros 
leaving either an equal number or fewer of differential equations as compared to the number 
of unknowns appearing as differentials. Therefore, while it did not appear to be so at first, 
a system in which some unknowns are present only as algebraic variables is actually a 
system of DAEs, and the number of boundary conditions required is fewer than first 
appeared. 

Systems that are obviously DAEs from the start also are not straightforward in 
determining the number of boundary conditions required. For instance, suppose the first 
order system of Eq. (3-33) has been manipulated to a form such that the number of rows of 
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A is equal to the rank of A. This case would appear to require as many boundary 
conditions as the rank of A. However, it is possible that the algebraic equations introduce 
an interdependence among the differential variables that reduces the rank of A. These 
examples are given simply to present the complexity that is involved in evaluating a system 
of DAEs in preparation for their solution. 

3.6 Resolving the Order of the System and Number 
of Boundary Conditions 

For the present case, that of the Euler equations resulting from the stationarity of the 
Hellinger-Reissner variational principle, the system of field equations and boundary 
conditions should be a consistent boundary value problem. However, at first inspection 
this is not obvious, and some manipulations of the system are necessary to arrive at a form 
that makes this apparent. 

As presented in Eq. (3-15), the principle contains 29 variables per layer (sixteen stress 
quantities, seven weighted displacements and three displacements at the bottom and top of 
the layer) for 29 N total. The condition that displacements be continuous between layers 
reduces this number by 3(N— 1), leaving 26A+3. There are also 26/V+3 terms in the integral 
composed of a varied quantity and its coefficient Of these coefficients 14 N are differential, 
composed of all IN of the equilibrium and IN of the compatibility terms. There are also IN 
boundary terms at each end. Substitution of the variations of the traction reciprocity 
conditions and continuity of displacements between layers reduces the number of variables 
by another 3 (AM) with the elimination of internal interface displacements. At the same 
time, the number of terms in the principle is reduced by another 3(AM) from the 
combination of the layer compatibility terms into the interface compatibility equations (Eqs. 
(3-24)). In this procedure, two differential compatibility terms combine to produce one 
differential compatibility equations at each interface (AM interfaces). Also, prescription of 
the x-direction traction on a surface combines two differential terms into one as is evident in 
Eq. (3-29). 

Thus, after invoking interfacial continuity and surface conditions the system has 23 N 
equations and unknowns with IN differential equations of equilibrium, either 6AM, 6 N or 
6AM-1 differential equations of compatibility (depending on the surface conditions), and the 
remaining equations are algebraic. In this form, difficulties arise in solving the system 
having to do with the number of differential equations and boundary conditions. Because 
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the number of differential compatibility equations are reduced without a like reduction in the 
number of differential equilibrium equations, it is possible to have a system with an odd 
number of differential equations. Application of the same boundary conditions at each end, 
giving an even number of boundary conditions, appears to be in conflict with this. This 
inconsistency was remedied by Pagano by using an involved interpretation of what he calls 
end conditions. He arrived at these end conditions by assuming the integrand of the 
functional vanishes at the ends. These conditions were then incorporated into his central 
difference approximation using three-point forward and backward differences at the ends. 

As an alternative to Pagano’ s interpretation, a simple manipulation of the differential 
equilibrium equations is used to reduce their number by the same amount as the differential 
compatibility. This was done by differentiating the first of Eqs. (3-23), eliminating o [\ ' ' , 
a^ +1) ' and < 7 1 < ‘ +1)/ by substituting from G[ k \ Gj* +1) and G, < * +1) and then substituting the 
result into G ( 6 k) . It is worth noting that the traction continuity equation that is differentiated 
here is the same one that combined two of the differential compatibility terms when its 
variation was substituted into the principle. It is also the only algebraic equation that 
contains only differential unknowns (i.e., they appear as differential terms in other 
equations). The new form of the equation is labeled & 6 k) and is given by 



(3-34) 


This equation is defined for k = 1 — > N - 1 . 

The ^-direction surface conditions require a similar treatment in order to avoid an odd 
number of differential equations. For u[ N) prescribed, a differential compatibility equation 
applies while for 1 prescribed, the applicable equation is algebraic: 

o[ 2 ] cos Yn - sin Yn = ?x2 (3-35) 

In a similar treatment to that used previously, Eq. (3-35) is differentiated, is 
eliminated by substituting from G ( 2 N) and the result is substituted into G ( 6 N) . This gives an 
algebraic form of G ( 6 N) , labeled G£ N) , and given by 
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or = -to ) ! <’ + [r; to - «»-. ) + to )']<’ - <’ + <’ + <’ + <’ 

+ 24 « - <’)+ 4[rito - *„->)- 2]<’ <3;r,6) 

— 2 JV-l) SeC Yn ^x 2 

Therefore, G 6 (W| = 0 is applicable when T^’ is prescribed while G^' = 0 is applicable 
when u 2 N) is prescribed. Following a similar procedure for layer 1 and G, 11 gives 


m x) 


= [rofe - Zo) - (Zo) 2 ]^!' + (4) 2 < - Oa” + *32 - < + < 

+4[2 + 7o(z, - Zoj\<$i + 2z;(-a‘‘ ) + <) - (z, - Zo)«c y 0 


f(D' 

‘■jti 


(3-37) 


3.7 Summary of Governing Equations 


The governing equations are summarized below with the differential equations 
underscored and the algebraic equations not underscored. 

3.7.1 Compatibility Equations 

xZ = 7% = 7® =x£=x%=xl = 0 k = l^N (3-38) 


^nr=0; i- 1.3,4, 5, 6 

^ < 2 > +z^ 2 ) +(z;) 2 ^2 ) = 0 

^2 , + Z^2 , =0 


Jt = 1 -» W - 1 (3-39) 


^ 1 1 ) +z^ 1 , ) =0 

*2 ) +z;*5? ,as 0 
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3.7.2 Equilibrium Equations 


d k) = G[ k) = G {k) = G\ k) = G {k) = 0 


G {k) = 0 
G <* + n = q 


with either, for u[ l) prescribed. 


G' 11 = 0 


or, for r' 1 ,’ prescribed. 


G 5 (1) = 0 


as well as either, for u[ N) prescribed. 


G ( 6 N) = 0 


or, for t^ 1 prescribed. 


G 6 W = 0 


k, = \—>N 
k = 1 N - \ 


3.7.3 Interlayer Traction Continuity Equations 


K-^;<) = K +u -^< +1) ) 

(< 7 ^ - z' k O ( £) = « +1) - zX +1> ) 
« 2 ’ - z&i) = K +1) - ^< +I) ) 


k = \^N-\ 


3.7.4 Surface Conditions 

For layer 1, 


Xi\ + z'oXu - “i (1) = 0 


or 

cr^siny, -cr' l' cos y 0 = 


(3-41) 

(3-42) 

(3-43a) 

(3-43b) 
(3 -44a) 
(3-44b) 

(3-45) 

(3-46a) 

(3-46b) 
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(3-47 a) 


or 

Oa sin Yo ~ a 4 i cos Yo = *5* (3-47b) 

Xm “ w, (1) = 0 (3-48a) 

or 

of/ sin y 0 - (Jj 1 ]' cos y 0 = (3-48b) 

and for layer N. 

= 0 (3-49a) 

or 

cr^° cos y N - crf^ sin y N = (3-49b) 

^ + v< N) =0 (3-50a) 

or 

= C (3-50W 

■/l' + *{ K ’ = 0 (3-5 la) 

or 

<cos y/v -<sin r , = f^ (3-51b) 


3.7.5 Boundary Conditions 

The boundary conditions are affected by the imposition of traction continuity and the 
surface conditions as well. In Eq. (3-15), there are 14JV boundary terms at the x = constant 
edges (C) of which IN apply for a particular set of conditions. Substitution of the 
variational form of the first of Eqs. (3-23) reduces the number of displacement prescription 
equations by N-l. Applying the unvaried form of the same equation reduces the number of 
independently specifiable stress variables in the traction prescription terms by N— 1 . In a 
similar manner, for each ^-direction surface traction prescribed, the number of boundary 
conditions is reduced by one, thus retaining a consistent number of differential equations 
and boundary conditions. 

A further point needs to be made regarding the development of the boundary 
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conditions. In order to allow the three tractions to be applied independently, restrictions 
must be placed on the geometry of the domain at the ends. Specifically, the slope of all 
layer interfaces are restricted to zero at the ends. The reasoning behind this restriction is 
illustrated in Fig. 3.5 where it can be seen that the ^-direction traction continuity equation 
imposes an interdependence between z x and T, if z[ * 0. 



Fig. 3.5 Enlargement of layer interface region at an end illustrating interdependence 

between traction B.C.s. 

The resulting boundary conditions are 


II 

?*• 

and 

u (k) = u (k) 

k = l->N 

(3-52a) 

or 





— * 

^tT 

ii 

* — . 
V 

and 

(J — Q-W 

U \2 U 12 

k = l -+N 

(3-52b) 

y<*) - y(k) 

and 

y(k) = y<*) 

k = 1 N 

(3-53a) 

or 





< = 

and 

<j (k) — 

U 62 U 62 

k = l—>N 

(3-53b) 


(z,., -z,)W«* u 




= (z>-z>-,W k ’ 

and 

+ (z».,-z,)W M> 


k = l N -l 

(3 -54a) 

W ik) = W ik) 



k = l^ N 


or 
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r <*+l) _ 0.(*+l> 


'51 


51 


and 


r <*) _ *<*> 


'53 


= < 7 . 


53 


With, for u[ l) prescribed, 


and for u[ N) prescribed, 


k = l->N-l 
k = \->N 


(3-54b) 


W w = w (l) 

(3- 55a) 

or 


C £ 

*b 

it 

(3-55b) 


(3- 56a) 

or 


A 

II 

Q* 

kjTW 

^5 

(3-56b) 


3.7.6 Continuity Between Longitudinal Segments 

In order to analyze a dropped-ply laminate, the formulation must also allow step 
changes in material properties in the longitudinal direction. To accomplish this, the domain 
is divided into segments, each having its own properties. The solutions for the different 
segments are joined through inter-segment continuity conditions. In deriving these 
conditions, restrictions on the geometry of the adjoining segments are imposed. 
Specifically, the layer boundaries must be continuous and the slope of the layer boundaries 
are required to be continuous. The first restriction is required to preserve consistency of 
definition of the weighted displacements between segments and to allow pointwise traction 
continuity. The second restriction is imposed because without it, Ci, a 3 and <7 5 would all 
have to be equal across both the layer and segment interfaces at the junction. The resulting 
conditions between segments designated a and b are 
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y(k) _ y{k) 


y(k) _ y(k) 
r a r b 

w (k) = w a) 

a ff a 


uw' k) 


u: 


<*) 


+ Z& 


(k) = U ( b k) + z[W b 


G (k) — / t <*> 
U 6\ a U 61 b 

.<*) _ ^k) 
62 b 


<nv* = <? 


'62 a 
(j(k) — (k) 

U 53 a U 53 b 

(j(k) _ q- ( k) 

U n a U 11 b 
u 12a u 126 


(*) 


— (fc+ 1) _ _(* + !) 


■(*) 


7 51a 


= a; 


516 


With, for Mj (1) prescribed, 

W a m = W b (U and 

and for u[ N) prescribed, 

W a (Af) = W fr (JV) and 


k = \ —> N 


k = l-> N-l 


<T (1) = <T (1) 
U 51a U 51fc 

„iN) (AO 

°52 a u 52b 


(3-57) 


(3-58) 


(3-59) 


(3-60) 


3.7.7 Accounting of Variables and Equations 

This completes the presentation of the equations associated with the model. The model 
is summarized in Table 3.1 by listing the unknowns and the equations used to solve for 
them. The number of first order ordinary differential equations, and the associated 
boundary conditions, are summarized in Table 3.2. Note that the order of the system is 
dependent on whether the .r-direction displacements on the top and bottom external surfaces 
are prescribed or not prescribed. 
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Table 3.1 Summary of the mathematical model (N = number of mathematical layers). 


Unknowns 

Type 

Number 

G\\ * ^12 * ^21 * ^22 * ^ 31 * ^32 * ^ 33 * ^34 » 

stress variables 

16 N 

^41 ’ ^42 * ^43 ’ ^ 51 ’ ^52 ’ ^53 ’ ^61 ’ ^62 



U,u,v,v,w,w,w 

weighted displacements 

IN 


Total 

23N 

Equations 

Type 

Number 

(3-38) 

layer compatibility 

6N 

(3-39) 

interface compatibility 

7(N-1) 

(3-40) 

surface compatibility 

4 

(3-41) 

layer equilibrium 

5N 

(3-42) 

interface equilibrium 

2(AM) 

(3-43), (3-44) 

surface equilibrium 

2 

(3-45) 

traction continuity 

3(AM) 

(3-46) through (3-51) 

surface conditions 

6 


Total 

23N 
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Table 3.2 Summary of first order differential equations and boundary conditions (BCs). 

N = number of mathematical layers. 


prescribed surface 
displacements; jc, < x < x 2 

Differential Equations 

BCs at X\ and x 2 


u[ N \x) 

Equations 

Number 

Equations 

Number 

no 

no 

(3-38) 

(3-39) 

(3-40) 

(3-41) 

(3-42) 

Total 

N 

5N-5 

4 

5N 

AM 

12N-2 

(3-52a) or (3-52b) 
(3-53a) or (3-53b) 
(3-54a) or (3- 54b) 

Total 

4^V 

4A^ 

4N-2 

12N-2 

yes 

no 

(3-38) 

N 

(3-52a) or (3-52b) 

4N 



(3-39) 

5N-5 

(3-53a) or (3-53b) 

4N 



(3-40) 

4 

(3-54a) or (3-54b) 

4N-2 



(3-41) 

5N 

(3-55a) or (3-55b) 

2 



(3-42) 

AM 





(3-43a) 

1 





(3 -46a) 

1 





Total 

12N 

Total 

\2N 

no 

yes 

(3-38) 

N 

(3-52a) or (3-52b) 

4N 



(3-39) 

5N-5 

(3-53a) or (3-53b) 

4N 



(3-40) 

4 

(3-54a) or (3-54b) 

4N-2 



(3-41) 

5N 

(3-56a) or (3-56b) 

2 



(3-42) 

AM 





(3 -44a) 

1 





(3 -49a) 

1 





Total 

\2N 

Total 

\2N 

yes 

yes 

(3-38) 

N 

(3-52a) or (3-52b) 

4N 



(3-39) 

5N-5 

(3-53a) or (3-53b) 

4N 



(3-40) 

4 

(3-54a) or (3-54b) 

4N-2 



(3-41) 

5N 

(3-55a) or (3-55b) 

2 



(3-42) 

AM 

(3-56a) or (3-56b) 

2 



(3-43a) & (3 -44a) 

2 





(3-46a) & (3 -49a) 

2 





Total 

12N+2 

Total 

12N+2 
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Chapter 4: Axisymmetric Model 


The next version of the model which will be developed is for the axisymmetric 
response of a shell of revolution. This is a reformulation of Pagano’s 45 theory for this class 
of problems with the incorporation of the changes applied in the generalized plane 
deformation model discussed in Chapter 3. In addition, a further modification was made to 
the development of the stress assumptions which reduces the number of dependent 
unknowns per layer by two from Pagano’s original derivation. 

The coordinates chosen for this derivation are slightly different from those used in most 
cylindrical coordinate systems. Specifically, the axial coordinate will be x rather than z, 
while the radial and circumferential coordinates will be the conventional r and 6, 
respectively. This change was made in order to have x as the independent variable as in the 
generalized plane deformation model and also to avoid confusion that may result from 
another z coordinate at a different orientation. 

In the present class of problems, the loading, geometric and material properties, and 
hence stresses and strains are independent of 6. The solution domain is thus reduced to the 
x-r plane. The stress field will be assumed explicitly in terms of r and implicitly in terms of 
x with independence of 6. Therefore, curvilinear coordinates (a,,a 2 ,£) used in Sec. 2.2.2 
are identified as (x, 6,r ) for this formulation. This model is applicable to the analysis of 
laminated cylinders with axially dropped plies as illustrated in Fig. 4.1. 

4.1 Displacement Field 

Unlike generalized plane deformation which allows for dependence of the 
displacements on the y coordinate, axisymmetry does not permit dependence of the 
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displacements on 6. In addition, there is no rigid body translation in the radial and 
circumferential directions. Therefore, neglecting rigid body displacements, the 
axisymmetric form of the displacement field is 

u(x,d,r) = U(x,r) 

v(jc, 6, r ) = V(x, r) (4- 1 ) 

w( x,6,r) = W(x,r) 



The engineering strain field for the displacements given by Eqs. (4- 1) has the following 
form: 


£ , = = u .* 

1 w W 

£ e “ v ,e ” 
r r r 

£ r = W, = 

v 1 V 

U=Vr— + -”e = V r — 
r r r 

Yxr = « , + w, = U', + W x 

=~ u ,e + v = v ,x 
r 


(4-2) 


4.2 Assumed Stress Field 

A generic layer for the axisymmetric model is shown in Fig. 4.2. The assumed stress 
field within such a layer for this model is taken in the form 
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(4-3) 


cr(x, r) = o u {x)f { j'\x,r), no sum on i 
where contracted notation is employed once again. 



Fig. 4.2 Schematic of a typical layer. 

As mentioned previously, the stress assumptions for this model were changed from 
those used by Pagano. These changes include those utilized in the generalized plane 
deformation model as well as an additional modification that reduces the number of 
dependent unknowns. Specifically, this additional change involves the assumption of the r- 
surface stresses <7*, and o x9 (<7„ cx 2 and <J 6 in contracted notation). In both the 
generalized plane deformation model and Pagano’s version of this model, the r-surface 
stresses were assumed to have a linear dependence on the thickness coordinate within each 
layer. This approach works well for the generalized plane deformation model giving the 
dependences of the other stresses shown in Eqs. (3-5). These expressions for the stresses 
have the characteristic of containing the minimum number of coefficient functions of x 
introduced through integration of the equilibrium equations. By contrast, if linear 
distributions of the r-surface stresses are substituted into the axisymmetric equations of 
equilibrium, given by 
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(4-4) 


~( r<7 tr) r + &XX.X = 0 

\(r 2 (J re ) r + o xe , x =0 

V) r +^,-- = o 

r ' r 

in which body force terms are neglected, the resulting expressions for the transverse 
stresses are 


1 2 , 

o n (x, r) = a x ( x ) - + b 2 ( x ) + c 3 ( x) r + ( x) r +e i (x)r 

r 

o & (*, r) = a 4 (x) \ + b A (x ) r + c 4 (x ) r 2 (4-5) 

r 

1 2 

O xr (*. r) = a 5 (x)- + b 5 (x)r + c 5 (*) r 
r 

Note that the expression for contains five terms rather than the four terms present in the 
generalized plane deformation model for o u (see a 3 in Eqs.(3-5)). 

Examination of the axisymmetric equations of equilibrium revealed that an alternative 
choice for the r-surface stresses would result in a form for a„ that contains only four 
terms. Division of the linear functions assumed for and <7 xe by r resulting in rational 
functions was found to achieve the desired result. This was not done with the third r- 
surface stress component, , because the 1/r term would produce a logarithmic term in 
Summarizing the general form for the assumed stress field in contracted notation: 

ct, (x, r) = c„ (x, r) = a, ( *)-J- + b } ( x) 

<T 2 (x,r) = <J 90 (x,r) = a 2 (x) + b 2 (x)r 

1 2 
oAx,r) = o n {x,r) = a i (x)- + b- i (x)+c 1> {.x)r + d i (x)r 

r i (4-6) 

<T 4 (x, r ) = O ffr (x,r) = a A (*)-j + b 4 (*)+ c 4 ( x) r 

cr 5 (jc, r) = <7 (x, r) = a 5 (at) ■ - + b 5 (x) + c 5 (x) r 

r 

<J 6 (x, r ) = <J xe (x, r) = a 6 (x)- + b 6 (x) 

r 

Once again, the final forms of the stress assumptions given by Eq. (4-3) are arrived at 
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through application of the conditions on the stresses at the surfaces of the layers given by 

&i (■** r \ ) = (x) ( 4 - 7 ) 

G i (x,r 2 ) = O i2 ( x) 

Also, conditions similar to those used in the derivation of the bubble functions for the 
generalized plane deformation model will be applied. 

The r-surface stresses are uniquely determined by the conditions of Eq. (4-7) and are 
given by 


a, =cr, 


a 2 = o 


21 


n(r 2 -r) 

r ( r 2~ r i\ 


+ CT, 


12 


r 2 ( r - r i) 
r(r 2 - /j ) 


+ a 


&(, ~ ^61 


v r 2 - n j 


r r\ij2~r) ' 


r-r. 


2 '1 > J 

\ 


22 


Ty-r, 


(4-8) 


+ o t 


62 


1 J 

/ V 2 (r-r 1 ) > 

r(h-r x )j 


\r{r 2 -r x ) 

The transverse shear stresses o 4 and <r 5 are developed through elimination of the 
coefficients b, and c„ i = 4, 5, by application of Eqs. (4-7) and subsequent use of the 
following expressions for a 4 and a 5 , which are chosen to scale and nondimensionalize the 
shape functions. 


a 4 = - 


2 2 


\2 ^43 


(r 2 -r x ) 

5 (r,-r, f 


(4-9) 


This leads to the following assumed forms for cr 4 and a 5 . 


&4 “ ^41 


. r 2 - y 


+ a, 


42 


<J 5 = Gi 


51 


^2 ^ 

\ r 2~ r i J 


+ <x 


52 


r-r 


V r 2 - r i J 


r-r, 


V r 2 “ r i y 


+ cr 


43 


+ a 


53 


-r 3 (r, +r 2 ) + r 2 (fj 2 + /jr 2 + r 2 2 ) - r t V 2 2 A 

r\r 2 ~r x f 

-(r 1 +r 2 )[r 2 -r(r,+r 2 ) + r,r 2 ] 
r(r 2 -r,) 2 


(4-10) 


The development of the assumption for (% is obtained in a different manner for the 
axisymmetric model than for the generalized plane deformation model. Instead of 
eliminating a couple of the coefficient functions through application of Eqs. (4-7), a desired 
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form for the stress assumption will be chosen and verified as having the form of Ch in Eq. 
(4-6) through the definition of all four of the coefficient functions. This change was made 
because the second bubble function obtained using the previous approach did not have all 
three of its zeros in the range r x < r < r 2 as is desirable to accentuate the difference in its 
contribution to the stress distribution from the first bubble function. 

The first two shape functions will simply be the familiar linear functions used in o 2 . 
The first of the two bubble functions is a scaled rational function arrived at by dividing the 
quadratic that passes through zero at r x and r 2 by r. The second bubble function is a scaled 
rational function composed of the cubic that has zeros at r x , (r x +r 2 )/ 2, and r 2 , divided by r. 
In each case, the scaling factor is chosen such that the peak value of the bubble function 
within the layer does not vary with r x and r 2 . This desired form for < 7 , results from the 
following definitions for the coefficient functions: 


_ r L r 2 (r t + r 2 ) -I- r 2 ) 2 

U 3 ~ X 2 °33 * 


( r 2 ~ r i) 


/ \3> w 34 

( h-r . ) 


- r 2 o 3x -r x o i2 (r, + r 2 ) 2 (fi +r 2 )(i i 2 +r 1 f 2 +r 2 2 ) 

(r 2 -r,) (r.-r,) 3 


'34 


C 3 “ 


_ ^32 z g 3 i (r^r 2 ) 3(r t +r 2 ) 


2 ^33 


( r 2“0 (^2 ^ 1 ) ( r 2 “' i ) 


x3 v 34 


d =^ r i + ^ G 

3 ( r 2 ~r x f 34 


and is given by 


(4-11) 


f . 


(T 3 - <T 3l 


r 2 -r 


+ a 


( 


r-r 


\ r 2^ r iJ 


+ <y 


+ a 


34 


~0j + r 2 )[r 2 - rjr\ + r 2 ) + r,r 2 ] 
r(r 2 -r x ) 2 

(j\ + h )[2r 3 - 3 r 2 (r t + r 2 ) + r(r L 2 + r,r 2 + r 2 2 ) - r,r 2 (r, + r 2 )] 

r(r 2 -/j ) 3 


(4-12) 


The shape functions of Eq. (4-3) are therefore given by 


A D = f (6) = *i(r 2 -r) 
* ^ r(r 2 -r.) 

f(i) _ f (6) _ r 2 (r-r,) 
r(r 2 -r,) 


(4-13) 
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^=f^/r=f ] = 


h~L 


-n 


f (2) _ f(J) _ f( 4 ) - f( 5 > = L — 

h ~ h ~ h h 


h ~' 1 


/f=/ 


(r 1 +r 2 )[r 2 -r(r 1 +r 2 ) + r,r 2 ] 


r(r 2 ~n) 


/i 3) = 


f(4) 

V3 


(r, + r 2 )[2r 3 - 3r 2 (r, + r 2 ) + r(r x 2 + r,r 2 + r 2 2 ) - r,r 2 (r, + r 2 )] 


r(r 2 -r,) 


-r 2 (r, + r 2 ) + r 2 (r 2 + r,r 2 + r 2 2 ) -r, 2 r 2 
r 2 (r 2 -r,) 2 


(4-13) 

(cont.) 


While this formulation produces an increase in the number of different shape functions 
from four in the generalized plane deformation model to seven here, it still represents a 
reduction from the nine in Pagano’s formulation. More significant is the reduction in the 
number of stress variables from seventeen in Pagano’s formulation to sixteen in the present 


one. 


4.3 Application of the Principle 

Because of the independence of the integrands of Eq. (2-7) from 6, integration with 
respect to 6 for any axisymmetric volume simply introduces a factor of 2n, which can be 
canceled, and reduces the volume and surface integrals to area and contour integrals in the 
x- r plane. Integration with respect to the thickness coordinate, r for the present model, 
again requires the application of Leibnitz’s theorem for terms involving derivatives in x and 

also the definition of weighted integrals for the displacements. 

The definitions of the weighted displacements for this model were chosen in a similar 
manner to those of the generalized plane deformation model in that they are based on the 
shape functions used in the stress assumptions. However, in this case, the shape functions 
are multiplied by r from the volume integration in cylindrical coordinates. Once again they 
are scaled by division by the layer thickness (r 2 - rj. Additionally, because of 
multiplication by r from the volume integration, they are also divided by either r„ r 2 or 
(r t + r 2 ) , based on which of these factors simplify the particular shape function. Finally, 
because of the changes in the stress assumptions, the same diacritical notation for two 
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displacement components does not necessarily mean that they are weighted by the same 
function. The weighted displacements are defined by 

(t7, U, V, V, W, W , W) = f r : (^U^U^V^V^W^W^W)^ (4- 14a) 

where the weighting functions are 




^2 — 
K = 


K = 


( r 2 - r \ ) 

lzJ± 

(h -'i ) 2 
_ r(r 2 - r ) 

r i( r 2~ r i) 2 

rjr-rj 

r 2 {r 2 -r l ) 2 


(4- 14b) 


(r 2 -r(? ! + r 2 ) + r,r 2 ) 
5 (r 2 -ii ) 3 


It is worth noting here that there is one fewer weighted ^-displacement unknown than is 
required by Pagano’s formulation. This is due to the changes in the stress assumptions 
which eliminate the need for an r 3 term in <%. 

Using the same notation for denoting layers and interfaces as was used in the 
generalized plane deformation model with the exception that r replaces z as the thickness 
coordinate, the variational principle for this approximation is 

[ 2 1 | [ 0 Uu + Xu r s °u k) - ( G J U + G 2 SU + G,SV + G 4 SV + G 5 SW + G 6 SW 

+G 7 5w) < * , ] + ([(r, 1 - t,,) 0 ’**}” +(t„ - T ei ) (n 5v; 1, +(T rl - f rl ) < "<5w 1 (1) ]r 0 secy 0 

+[(*,2 - T x2 ) <iV) 5Mf ) + (t 82 - r e2 ) (N) 8v[ N] '+(T r2 - i r2 ) (N) sec Yn) c ^ (4 ’* 

-([(«! — -t- (vj -wJ^Wjibsccyo 

+[(«2 +( v 2 - v 2 ) iN) 8r t B2 ) +(w 2 - w 2 ) (N) ]r w sec 7 n 
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£[« - ogv;- <*" + <*" v/)&4*’ +«’ - <’<- ’ + <' v,')^ 

*=1 

(oS$> - &£r;- <" + <VK]i V* 


+ Z{([ r *-i( CT '> + r *( CT « - a»f« w +U* - 

*i(ob - <r H f 5 V<‘> + U«« - 6 st fSW'» + *(«* - a^f «V?“’ 

+(<U + r„)(a s , - a^fWfa - r *-.)) c , - (M* 7 - 


(4-15) 

(cont.) 


r k [u - 0^ Serf? + r*_,(v - v) ( ’&C + r*(v - S( + r*_,(w - w) Sag 

4*-*r**^t*-*i >5<y » ~ r *-i)) c j = ° 


It is worth noting that, except for the appearance of the radial location of layer interfaces 
multiplied by some of the terms, Eq. (4-15) is the same as Eq. (3-15) for the generalized 
plane deformation model. Therefore, the terms will not be described again here and the 
quantities appearing in the principle will be presented. 

The surface displacements are contained in Hu which are 


jC = r *-i r *-X‘ > 

= ~ r k-A k) 


Hn=- r k < u V 

< 2 ’ = r k W 2 k) 

Ha2 = 


< = + £,<') < = <iK’ -’■>?’) 


!«(*) — y y f 

A*61 r k-\'k-l v \ 


Ha =-h< v i k) 


(4-16) 


The Xu are defined similarly to the generalized plane deformation case with the exception 
that the out-of-plane prescribed strain terms do not appear in the present model. 


y(*) _ «(*) _ stW 

Xu ” r lu °ijKJ U jK 


(4-17) 


The weighted compliances are defined by 


?(*) 

^ijKJ 


r ‘ Sij fx'ff'rdr 


(4-18) 
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and the strain weighted displacement quantities are 


<’ = -h-,)u m ' +r t -A<- 2r U)u m -'i-i'-.'t/'” 

n\\' = -i(r, - r, t7“’ + r,(2r,'-r, 

lit’ =('■» -[(<•* 

'll;' = (2r„, -r,)iv» + +[(r, - r,.,)7r,_,]lV" 

< = -r l _,r*' + (r 1 . 1 -2r,)t»'‘*’ -[(r, -r,.,) 2 /,;]** 

C = -('i-, + r k )W<" + (r,., + rjH' 11 ' + < r ‘- + r »> ( r « ~ r » d) l v'>' 

< = -(>!-, +-i)r < » 1 -(r._, + r k )W m + (r *- + <i ^‘ - W'*’ 

7l“> - (3r M - 2r,)F'*’ + r t V< 1 ' < = + (2r,_, - 3r t )V" 

+ 2r,)v<*> +(2r t _, + rJV”' 

< = r i-i{ r k - ’■k-,W U, '+( 2'i-i - -iF*’ + rfi 1 

A’-k-,«-X-,)+v;_,]W' i> - v ;w'‘ 

< = r,(r, -r M )W*>'-r M lT'*' + (r,. 1 -2r,)0" 

+'i_,r t '_,W“ > + [r,(3r,'- r,'_,) - r,_ ; i 1 ']VV , ‘ 

< =k., +r,)[(r, +0'*' H'.if 1 *’ 

'i-i’iti'i* + Vi r i + r > 2 -i ) ~ 'i’i-i ( r i 2 + 3r »-' , i + 2r » 2 -' ) 


/(*) 


/(*) 


>(*) 


;(k) 


y(k) 


j(k) 


j(k) 


■;v (k) 

ik) 


C = U'i -r 1 . 1 )V«’' + r,.,( r ;- 2r/.,)V“> - r^V 
< = r *( r * - r,- l )V'* l, +™'-,V , ‘’ + r,(2r;-r,'.,)V' 

The equilibrium equations are 


( 4 - 19 ) 
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Gl k) = 


G ( 2 k) 


a k) = 


G<*> = 


Gf 


G (k) = 


G|*’ = 


.(*) 

'52 


T (k) 

M3 


+0*-l+»i) O » > 

+(2r t -r t _,)(y 5 

"fo-i + r k )<% 

f k -i( r k ~ r*_, X'+fcr* - 3r*_, )<> + r t _X + (r 4 _, + 2r k )& t 
+r k £ l (o%'-og) 

r k {r k - r t _i)og )/ -r t < > + (3r t - 2r t _ 1 )o£ ) -(2r t _, + r t )<r‘ 

+'i-i^« > “ < 2 ) 

r k -i{r k -»i_ 1 X )/ -(r* -r k _,)a^ + (r k -2 + r^og 

+(r*_i + »i X°» + °») + r *-i r *-i K’ “ <’) “ r *-i( r *-i + r * )°» 

n(r t - r t-i) <7 52 >,_ ( r * - r k _ x )o% - rf +{2r k - r k _ : )o£ 

~{r k . x + r k )(o l £ -<’) + r t r;(<> - a£’ ) + + r t )< 3 ’ 

(r , 2 - 'il,K 1 '+ ( ' i ~' i " -('i ^ 1 - 'i-,< - «i< + 

r it-i r jc 

- ( —-- r ^[(r t 2 - £,)< + (r k + 4r k r t _ 1 + i^X] 

r k -\ r k 

, ^<-i(3^-, + ^)~ + 3r t ) 


.(*> 

'43 


(4-20) 


■a (k) 

°S3 


Due to the changes made in the stress assumptions from those that Pagano used, the 
present axisymmetric model has the same form as the generalized plane deformation. That 
is, there are an equal number of terms in the variational statement and the same terms 
contain the same differentiated variables. Therefore, the manipulations of Sec. 3.4 are the 
same as those used in the present model and will not be reiterated here. 

The equilibrium equations resulting from the manipulations are as follows: 
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& 6 k) = ~r k -\{ r kf + [ r k ( r k - r k-i ) r " + r *-i (O^K ~ ( r * " r *-> ) G n 

-*< + (2r k - r w K’ " (*-. + OK " O 

+r*'[2r t < ) -(3r t -r k _,)&£ + %-i + 0°£ > ] 

+ r z M {[ r *-'( r *) 2 

\ r k+\ r k) 1 

('if -'•»)<’ -fo.. -2>iK*" 

-r,<" -k + r„, )«■’ + <") 

+tf(r w -3r,)<" + 2r,<» + 2(r ll +-;.,)<"]} 

which is defined for k = 1 — * /V - 1 . 

G 5 (1) = [rjr.(r, ~ r 0 )- r,(r 0 f]< + -(r, - r 0 K» + (r, - 2r 0 )< 

+^0 <J 32 + ( r 0 + X°» + <T w) + r o\r'o r O ( r i ~ r o) + \- r i K 
-2tf[r 0 Oja + (r 0 + r,)^] - r 0 (r, - r 0 )sec y 0 r^' 


C 6 (W) =-r A ,. 1 fe) 2 <’+[r^fe -r Ar . I ) + r Ar .,(r j ;) 2 ]< W) - (r # - r*_, Jo**’ 

~ r N G l\ ) + {2 r N ~ r N-l) (J y2 ) ~( r N-l +r A/)( a 33 ) _ °34 ’) 

+2r;[r„<° + (r w _, + r J ,)o^ , ] + ^[^(r w -3r N ]cr^ ! 

+0 r * — fy.Jsec Y fj t x2 


4.4 Summary of Governing Equations 

The governing equations for the axisymmetric model are summarized below. Once 
again, the differential equations are underscored and the algebraic equations are not. 

4.4.1 Compatibility Equations 

A? = x£ = Xyh = = X% = X% = 0 k = 1 -> N (4-24) 
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XI 2 + Xn +l) = 0; / = 1,3, 4, 5, 6 

z£' + r;ti»+{rffz*' =0 


k = l^N-l 


Xu +'0X5" +{ r of X* =0 

x£+tol = 0 
zS'+r'stiZ' + tifzg' ^O 
z£ ) +'izS ) = o 


4.4.2 Equilibrium Equations 

Gl k) = G'*’ = G<*> = G<*> = Gf = 0 


G‘*> = 0 
Gj* +1) = 0 


with either, for u[ ’ prescribed. 


Gj n = 0 


or, for t”/ prescribed. 


G 5 U) = 0 


as well as either, for u ( 2 N) prescribed. 


G[ N) = 0 


or, for r { X 2 prescribed. 


G 6 (A0 = 0 


k = l-*N 
k = 1 -> N-l 


(4-25) 

(4-26) 

(4-27) 

(4-28) 

(4- 29a) 

(4-29b) 
(4- 30a) 
(4-30b) 
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4.4.3 Interlayer Traction Continuity Equations 


«-«) =«"-«") 

« - '■;<) = «*" - 


k = l^N-l 


4.4.4 Surface Conditions 


For layer 1, 


or 

On’ sin 7 0 - cos y 0 = 


or 

< ) siny 0 -< ) cosy 0 = T; l 1 ) 

=0 

or 

^si’ sin To ~ <T 3i > cos y 0 = rjj* 


and for layer N. 


x^+r'xi r+«T' = o 


(W) , -<A0 _ , 


or 


<*s? cos y N - off sin y* = t<*> 


^+vr=o 

or 

<C cosy* - sin y„ = 


(4-31) 

(4- 32a) 

(4-32b) 

(4-33a) 

(4- 3 3b) 
(4- 34a) 

(4- 34b) 

(4-35a) 

(4- 35b) 
(4- 36a) 

(4- 36b) 
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4.4.5 Boundary Conditions 

The boundary conditions are derived in the same manner as for the generalized plane 
deformation model. The conditions on the stress variables and their variations required by 
the jc-direction traction continuity equation reduce the number of independently specifiable 
boundary conditions by 2AM. In addition, the interface between layers is restricted to 
having r - 0 at the ends of the domain as before in order to eliminate interdependence 
among the three prescribed end tractions. 

The boundary conditions are 


U lk) = u {k) 

and 

q< 

11 

«2« 

n- 

k = \^>N 

(4-38a) 

or 





fT (k) ~ & lk) 

Oil U 11 

and 

rr (k) - & k) 
°12 V 12 

Jt = l-> W 

(4-38b) 

y (k) = y {k) 

and 

y(k) _ y( *) 

k = \^N 

(4-39a) 

or 





II 

A 

25 

and 

n (k) _ gtt) 
°62 °62 

k = l-*N 

(4-39b) 

\r k -r k . ,)W<‘> + 

(v,-r.)<r- > 




and 



fc = l->N-l 

(4-40a) 

II 



Jt = 1 — > AT 


or 





' (ik+D _ Pr ik + l) 
°51 “ U 51 

and 

O53 - o 53 



jfc = l-> N-l 
k = l-*N 

(4-40b) 


With, for kJ” prescribed, 
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(4-4 la) 


and for u ( 2 N) prescribed, 


W (1) = W (1) 

or 


fT <U = <T ( " 

U 51 U 51 


^ (JV > = 

or 

< 7 ,N) = a {N) 
U 52 u 52 


(4-4 lb) 

(4-42a) 

(4-42b) 


4.4.6 Continuity Between Longitudinal Segments 

The conditions between segments of the domain also parallel those of the generalized 
plane deformation model. Once again, the layer boundaries and their slopes are required to 
be continuous. The conditions between segments a and b are 


y(*) 

a 

-X 

iS? 

n 


y(k) 

a 

= Yt k) 



= 


ul k) 

+ C 1 WJ k > 

= u^ + cx k> 

ul k) 

+ rX ] = 

U + r/Wf’ 


= G ik) 
U 61 b 


a (k) 

v 62 a 

= <7 {k) 

u 62b 


(T '(*) 

— 

U 53 b 


U \\ a 

= a (k) 
u n b 


>j(*) 

U 12 a 

= <; 


( r k + 1 


+ ( r * 


k = \^N 


(4-43) 




(*) 


U 51 a ~ °51 b 


& = 1 — > N - 1 


(4-44) 


With, for m, < 1) prescribed, 
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and 


( 4 - 45 ) 


w a n) = w b (l) 


(T* 1 * = CT*** 
°51 a U 51b 


and for u ( 2 N) prescribed, 

W< N) = Wl N) and < = <C 

This completes the presentation of the axisymmetric model. 


( 4 - 46 ) 
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Chapter 5: Numerical Solution and 
Verification of Accuracy 


Application of the models presented in the previous two chapters to laminates requires 
the definition of the geometric and material properties and the solution of a two-point 
boundary value problem. This chapter presents the methods that were chosen to accomplish 
these tasks as well as some that were found to be less suitable. In addition, the limits of the 
applicability and the accuracy of these structural models are examined. 

5.1 Modeling of Dropped-Ply Laminates 

The cross section of a dropped-ply laminate is modeled as shown in Fig. 5.1 consistent 
with the developments presented in Chapters 3 and 4. Because of the relatively large 
number of dependent unknowns per mathematically modeled layer, efficient 
implementation requires minimizing the number of layers used. This is accomplished by 
including multiple plies (or sublaminates) in each mathematical layer and by using thicker 
layers to model larger groups of plies located away from the stress concentration. 

Therefore, a method of determining sublaminate material properties consistent with the 
assumed stress fields is required. 

Another feature of dropped-ply laminates that must be accounted for is the material 
discontinuity at the end of the terminated plies. As mentioned previously, this is modeled 
by longitudinal segmentation of the model with different material properties in adjacent 
segments. The continuity conditions developed in Secs. 3.7.6 and 4.4.6 are then applied 
between these segments. 
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ZOT r 




Fig. 5.1 Modeling of a dropped-ply laminate. 

A final note that should be made regarding the model of the dropped-ply laminate 
shown in Fig. 5. 1 involves the layers used to model the terminated plies. Because these 
models are developed for bodies with continuous mathematical layers, those layers that 
constitute the dropped plies cannot simply be terminated with the plies. Therefore, they are 
continued through the triangularly shaped resin region and also through the thin portion of 
the laminate as very thin layers between the lower and upper continuous sublaminates. 


5.2 Incorporation of Material Properties 

The material properties of the mathematical layers are introduced into the models 
through the weighted integrals of the compliance coefficients S i]KJ defined in Eqs. (3- 1 8) 
and (4-18). For layers composed of more than one ply, these integrals can basically be 
determined in two different ways. The first is to integrate piecewise within a layer such that 
each individual integral is bounded by ply interfaces, and assume that the compliance 
coefficients are constant within each ply that comprises the layer. This allows the 
compliances, 5 1; , to be moved outside of the individual integrals that sum to the layer value 
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of s'£,. 

The other option is to determine some approximation to the effective compliances for 
the layer and to then carry out the integration in one step through the entire layer. This 
approach allows for varying the complexity of the approximation, and therefore, the 
accuracy attained. The simplest form of this type of approach would be to determine the 
compliances for an “equivalent” homogeneous layer through some form of approximation. 
This homogenization approach was chosen to compare to the piecewise integration of the 
weighted integrals of the compliances. 


5.2.1 Piecewise Integration 

Consider first the option of breaking the integral up into a sum of ply-by-ply integrals. 
This appears to be a desirable option as compared to homogenization because it accounts 
for the location of the individual plies through the thickness of the layer due to the 
multiplication by the shape functions. For instance, since the first two shape functions for 
o„ have their peak values at the layer surfaces, longitudinally stiffer plies being located 
there would be reflected in less compliant Sj£j terms that represent bending. Using this 
approach, the integral of Eq. (3-18) would be broken up as follows: 


r (A) 
°ijKJ 


= isj T f/ff'dz 

l-i 


(5-1) 


where M is the number of plies within the layer, 5‘° are the compliances for layer /, and 


l 

a i =vi+X*« 

a-\ 

with a M = z k . 

5.2.2 Determination of Homogenized Compliance Coefficients 

A simple homogenization scheme was chosen to compare to the piecewise integration 
approach. The method is presented in Appendix A and is based on a combination of the 
Voigt and Reuss approximations. 47 The Voigt approximation is applied to the in-plane 
strains and serves to impose interply continuity of the in-plane displacements. The Reuss 
approximation is applied to the out-of-plane stresses and serves to impose interply traction 

reciprocity. 


(5-2) 


/ = 1 — > Af 


Numerical Solution and Verification of Accuracy 


69 



The result of this procedure is a set of effective compliance coefficients for a 
homogeneous layer that is, in a strain energy sense, “equivalent” to the sublaminate. 
Because these compliance coefficients are constant throughout the layer, the weighted 
integrals of the compliance coefficients for the generalized plane deformation model can be 
written as 



with a similar result for the axisymmetric model. 

5.2.3 Comparison of Approaches 

A comparison between these two approaches can be made by a simple analysis of a 
uniform laminate loaded in uniaxial tension. A uniform thickness [±45/0/90] s quasi- 
isotropic laminate under uniform extension (£ n =0.1%), with transverse distortion 
restrained (£ =0.0), and unrestrained to z-direction contraction, was analyzed. The 
thickness of the laminate is 2.54 mm and the material properties of the orthotropic layer are 

E n = 128 GPa (1 8. 5 Msi), E 22 = E 33 = 1 1. 3 GPa (1. 64 Msi) 

G 12 =G B =6.0GPa (0.87 Msi), G 23 = 3.38 GPa (0.49 Msi) 
v n = v i 3 = 0-3, v 23 = 0.35 

This problem was analyzed with the generalized plane deformation model using both 
approaches to determining the weighted integrals of the compliances. The results are 
presented in Table 5.1 along with results from analysis of the same problem using classical 
lamination theory. 


Table 5.1 Comparison of methods of determining layer compliances. 



Piecewise 

Homogenization 

CLT 


Integration 

Scheme 


46.8 kN/m 

143.1 kN/m 

143.1 kN/m 

N y 

8.76 kN/m 

43.2 kN/m 

43.2 kN/m 


These results illustrate a clear discrepancy between the two approaches in analyzing this 
simple problem. This difference, as well as the similarity of the results for the 
homogenization routine and classical lamination theory, can be traced to the way in which 
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these approaches model the interaction between layers in the sublaminate. 

The homogenization scheme assumes uniform in-plane strains throughout in arriving at 
the average laminate properties. CLT, in general, does allow for non-uniform in-plane 
strains through the thickness. However, for this case of no bending, they are constant 
throughout the laminate and therefore the response would be expected to be very close to 
the homogenization results. 

The use of piecewise integration of the weighted integrals of the compliance 
coefficients, on the other hand, makes no assumption regarding the in-plane strains 
directly. Instead, this approach assumes that the stresses have a continuous distribution 
through the thickness represented by the assumed stress field used in the theory. The in- 
plane strains are related to the stresses through the constitutive relation for each ply and are 
therefore, in general, discontinuous between plies as illustrated in Fig. 5.2. Thus, the 
piecewise integration method is more closely analogous to a Reuss-type assumption for the 
in-plane stresses, which for a bonded laminate, is clearly a poor approximation. The 
homogenization approach was therefore chosen for determining the layer material 
properties for the remainder of the present analysis. 


z 


constitutive 


relations 



Fig. 5.2 Strain distribution associated with a uniform axial stress distribution. 


5.2.4 Material Properties for Curved Sublaminate Layers 

Since the theories presented in this research are able to analyze laminates with curved 
plies, the material properties for layers of curved sublaminates must be determined. The 
method used is very straightforward and is described here for the sake of completeness. 

The procedure determines the sublaminate compliance coefficients through the 
homogenization procedure as if it were flat, then rotates these compliances about either the 
y or 6 coordinate, depending on the model. Referring to Fig. 3.4, the angle of rotation 
used is an average of the angle of the lower and upper surfaces of the layer. 
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-(ft i + y* ) / 2 , which will be referred to as <p. The present application of these models 
does not account for the effect of variation in thickness of a ply on the material properties 
through changing fiber volume fraction. 



Fig. 5.3 Rotation from ply x, z' coordinates to structural x, y, z coordinates. 

Fig. 5.3 shows the structural coordinates x, >\ z and the rotation (p about y to arrive at 
the ply coordinates jc\ y\ z. For such a rotation, the stresses and strains transform 
according to the tensor transformation relations: 

o' = a^a^Ou 

and (5-4) 

e'j = 

where a is the direction cosine matrix given by 


cos (p 0 -sinp 
0 1 0 
simp 0 cos<p 


(5-5) 


Taking into account the symmetry of the stress and strain tensors, switching to engineering 
strains, and applying contracted notation leads to the transformation relations given by 
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(5-6) 


o'=T a o 

and 

e'=T £ e 


where 


T ff = 


and 


T = 


cos 2 <p 

0 

sin 2 (p 

0 

-2sin<pcos<p 

0 

0 

1 

0 

0 

0 

0 

sin 2 <p 

0 

cos 2 (p 

0 

2 sin (p cos <p 

0 

0 

0 

0 

CO s<p 

0 

sinqp 

sin (p cos (p 

0 

-sin<pcos<p 

0 

(cos 2 <p - sin 2 (p} 

0 

0 

0 

0 

-sin^p 

0 

cos<p 

cos 2 (p 

0 

sin 2 <j o 

0 

-sin (pcos<p 

0 

0 

1 

0 

0 

0 

0 

sin 2 (p 

0 

cos 2 (p 

0 

sin<pcos<p 

0 

0 

0 

0 

cos<p 

0 

sin i 

2 sin (p cos (p 

0 

-2 sin (p cos (p 

0 

(cos 2 (p -sin 2 (p ) 

0 

0 

0 

0 

-sin q> 

0 

cos 


(5-7a) 


(5-7b) 


Therefore, starting with the Hooke’s law relation in ply coordinates, 

e '= S'ct ' (5-8) 

substitution of Eqs. (5-6) leads to the following expression for the compliances in the 
structural coordinate system 

S = T; 1 S'T a (5-9) 


in which T~ (<p)=T e (-<p). 


5.3 Two-Point Boundary Value Problem 

Application of these models requires the solution of a two-point boundary value 
problem composed of a system of first order linear DAEs and their associated boundary 
conditions. This section presents the problem in general form and the solution method 
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chosen for the present work. 

In matrix form, the most general form of the system of DAEs to be solved can be 
written as 

A(jt)y'Oc) + BU)yU) = f(x) (5-10) 

where A(x) is singular due to the presence of the algebraic equations for which the 
associated row in A(x) is all zeros. The coefficient matrices A(x) and B(x) contain the 
geometric and material properties of the particular problem being analyzed. Due to the 
curved geometry at the ply drop-off, the system that the current work is concerned with has 
variable coefficient matrices as in Eq. (5-10). The solution vector y for an W-layer problem 
has the form 

y r =[y"’ r .y < 2 >r .— y w>T ] <5 - n) 

where each y ik) contains 23 unknowns composed of the 16 stress variables and 7 weighted 
displacements, 


v (k)T = fc 


( 5 - 12 ) 


It is worth noting that the A(x) and B(x) matrices have a block diagonal structure. This 
is due to the fact that each of the equations contain unknowns for, at most, two adjacent 
layers. While the present work does not take advantage of this fact, it could prove useful 
for improving the efficiency of the solution of the models. 


5.3.1 Solution Method 

The solution method chosen for this system is a one-step finite difference 
approximation. One-step differences were chosen over two-step schemes such as the 
central, forward, and backward differences used by Pagano 48 - 49 because of their flexibility 
in handling non-uniform mesh spacing as well as their lack of the need for special treatment 

at the boundaries for first order systems. 

The application of a finite element approximation in x was also considered. This 
approach was developed for Pagano’ s original laminate theory 44 by Sandhu et al. 50 This 
theory applies to flat laminates and results in a system of partial differential equations in the 
plane of the laminate. The procedure Sandhu et al. followed was to develop a self-adjoint 
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version of Pagano’s governing equations, which were not self-adjoint as published, from 
which they derived a new variational formulation. Using the Euler equations from this new 
variational formulation, the layer stress resultants were then solved for in terms of the 
displacement variables and the interfacial tractions through the constitutive/compatibility 
equations, thus guaranteeing their satisfaction. Finally, the displacement variables and 
interfacial tractions for all of the layers were approximated in terms of nodal values and 
interpolation functions in the 2-D domain of the laminate in an application of the finite 
element method. However, this solution method was very involved even for the flat 
laminate, constant coefficient case treated by Sandhu et al. Verification of self-adjointness 
of the equations for the generalized plane deformation and axisymmetric models is 
complicated somewhat by the fact that they have variable coefficients in contrast to the flat 
laminate formulation of Sandhu. If the governing equations are determined to be self- 
adjoint, they still are not formulated in terms of stress resultants, and because the layer 
interfaces are not required to be straight and parallel to the x-axis in general, five 
components of stress can contribute to the surface tractions as opposed to just three for the 
flat laminate case. This would make selection of interpolation functions for the interface 
stress variables such that they satisfy interlayer traction continuity much more difficult than 
in the flat laminate case. It is not clear what changes to Sandhu’s approach would resolve 
these issues and therefore, such a finite element approach was not considered practical at 
this stage of the development of the present models. 

In applying one-step finite differences to a first order system of ordinary differential 
equations, the approximate solution is represented by values of the dependent variables at 
discrete mesh points throughout the domain. Consider a finite difference mesh of <p points 
distributed along the x-axis. The location of each mesh point isx„, a = 1,2,...,0, and the 
solution vector at x a is y„. Solution by finite differences simply replaces the differential 
equations by a set of algebraic equations in terms of the mesh values of the dependent 
variables. One-step finite differences are a sub-class for which each algebraic equation 
involves the values of the dependent variables only at two adjacent mesh points, i.e. one- 
step. Examples of one-step difference schemes are the trapezoidal and midpoint schemes. 

In matrix form, the finite difference equations for a one-step scheme approximation to a 
system of p first order ODEs can be written as 

Q«y a +R a + iy B+ i =s a a = 1,2,.. .,0-1 (5-13) 

where Q and R are p x p matrices and s a is a p x 1 vector, all of which are determined for 
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the interval from x a to x a+l . This, along with the p boundary conditions at x x and x ¥ 

i>,yi + = e (5 ' 14) 

where D, and D, are p x p matrices and e is a p x 1 vector, gives p<p equations for the p<t> 
unknown values of the dependent variables at the mesh points. The matrix equation used to 

solve for the approximation is 


Qi 

r 2 


’ yi ' 


S, ' 


Q2 ^3 


y 2 


s 2 


*. 



► ~ < 

; . 


Q*-i 

R * 

y*-. 


Vi 

_ D > 


D *. 

.y# . 


e 


(5-15) 


In contrast, solution of the same system by two-step finite differences, which are 
usually applied to second order differential equations, is less definite. Application of central 
differences to a system of first order differential equations, as done by Pagano, results in 

the matrix equations 


P«-iya-l + Qaya + R a + iy« + l _S a 


a = 2, 3,.. .,0-1 


(5-16) 


which apply at each internal mesh point for p(^-2) equations. Boundary conditions such as 
those of Eq. (5- 14) provide another p equations leaving a shortage of p equations for a 
complete system. As mentioned previously, Pagano used a combination of two-step 
central, forward, and backward differences with the forward and backward differences 
applied at the first and last mesh points respectively. For a system of ODEs, this approach 
would apparently leave no room for incorporation of the boundary conditions. However, 
for a system of DAEs, manipulations can be made that, in a sense, incorporate the 
boundary conditions into the conditions at the end mesh points. Pagano accomplished this 
by assuming that the integrand of the variational principle vanishes at the ends of the 
domain and derived “end conditions” which were a combination of the boundary conditions 
and some of the field equations. Experimentation with this approach presented by Pagano 
led to difficulties due to a step change in the governing equations between the end mesh 
points (the end conditions) and the interior points adjacent to them (field equations) for 
certain displacement prescribed boundary conditions. Having resolved the number of 
differential equations and boundary conditions as described in Sec. 3.6, the decision was 
made to switch to a more straightforward one-step difference scheme. 
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Solution of a system of linear DAEs by one-step finite differences is similar to the 
handling of ordinary differential equations. In fact, the differential equations of the system 
are replaced by difference equations in the same manner as for the ODEs. However, for a 
system of DAEs, there are fewer differential equations than there are unknowns. The 
additional equations required are simply the algebraic equations evaluated at each of the 
mesh points, including the end points. 

For the two structural models presented here, there are 23N equations and unknowns. 
Depending upon whether the x-direction surface tractions or displacements are prescribed, 
there are either 12/V-2, 12 A, or 12/V+2 differential equations, the remaining algebraic (see 
Table 3.2). Correspondingly, there are either 6N-1, 6 N, or 6N+1 boundary conditions at 
each end of the domain. For problems that require segmentation because of discontinuous 
material properties, two mesh points are placed at the junction between adjacent segments, 
one associated with each segment. The additional 23N equations needed for these 
additional 23 N values of the dependent unknowns are a combination of the segment 
continuity conditions of Eqs. (3-57) through (3-60) for the generalized plane deformation 
model or (4-43) through (4-46) for the axisymmetric model, and the algebraic equations 
applied at the additional mesh point 

5.3.2 Failure of the Trapezoidal Finite Difference Scheme 

The trapezoidal finite difference scheme was attempted first on the DAE system from 
the generalized plane deformation model. This scheme gets its name from the fact that it can 
be derived by approximately integrating the differential equations from x a to x a+ i using a 
single step of the trapezoidal rule for numerical integration. It has an accuracy of order 2, 
meaning the error is 0(h 2 ), h being the mesh spacing. However, as will be demonstrated, 
this scheme lacks the accuracy necessary to solve problems that contain rigid body 
displacements (or problems with large domains for which local displacements may be large 
relative to the local strains). 

The implementation of the trapezoidal scheme used here is for a first order system with 
non-constant coefficients. For the system of Eq. (5-10), the finite difference equations for 
the interval from x a to x a+l are derived by approximating the integral of the system from x a 
to x a+ , using a single trapezoid. Starting with the integration, 

['"‘Ay 'dx+ ['"' Bydx= [°"tdx (5-17) 

J *a Jx a Jx a 
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which can be rewritten as 


P* 1 [(Ay)' - A'y]dx + By dx = \ f dx (5-18) 

Jx a L Jx « x « 

giving 

[Ay]^*' + (B - A')y dx = J°' f dx (5-19) 


Approximating the integrals by a trapezoid leads to the finite difference formula 


[(b„-a;) a„ 

WJ _L 

(®a+l ^a+l) j ^a+1 

fa+1 

2 (*«+!-*«). 

y« ^ 

2 (**a+l ^a) J 

J a+1 ^ 


which is used in the following example. 

A simple problem used to test the trapezoidal scheme is shown in Fig. 5.4. The 
problem consists of a rectangular plate of uniform isotropic material subjected to a rigid 
body translation. While this choice of a rigid body translation may seem peculiar, the 
reasoning behind it is to simulate the situation a large distance from a fixed point where 
displacements can be large relative to the strains. The plate is 1 cm thick and 10 cm long 
and is made of aluminum alloy with E = 68.95 GPa and v = 0.31. The rigid body 
displacement is a uniaxial translation in the z direction of w = 0.5 cm and is imposed at the 
ends. The plate is modeled by two layers separated by a curved interface, which is 
antisymmetric with respect to x = 5, as shown in Fig. 5.4. As required by the model in 
Chapter 3, the slopes of the surfaces and internal layer interface are zero at the ends. 

The results of this test of the trapezoidal scheme as applied to the generalized plane 
deformation model is shown in Figs. S.5-5.7 for three levels of mesh refinement. Fig. 5.5 
shows the axial stress o a along the internal interface at both sides of that interface for a 
mesh composed of 40 equally spaces points. It is clear that these results are unacceptable 
considering this problem should produce no stress while the erroneous values produced are 
a significant percentage of the tensile strength of this material, which is 90 MPa. Fig. 5.6 
shows results for the same problem with an increase in the number of intervals to 80. The 
Oih 1 ) accuracy of this scheme is evident in the reduction of the value of the erroneous 
stresses by a factor of four resulting from a halving of the mesh spacing. 
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Fig. 5.4 Rigid body displacement test problem. 
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Fig. 5.5 Test problem results for w = 0.5 cm translation, solution by trapezoidal scheme 

solution with 40 intervals. 



Numerical Solution and Verification of Accuracy 


79 




40 


30H 


top of layer 1 

— — - bottom of layer 2 



204 



-10H 



T 

4 


x 

6 

x (cm) 


T I 

8 10 


Fig. 5.6 Test problem results for w = 0.5 cm translation, solution by trapezoidal scheme 

solution with 80 intervals. 

It appears from these results that simply increasing the number of mesh points will 
eventually produce acceptable results. However, increasing the number of mesh points to 
320 produces the results shown in Fig. 5.7. While the magnitude of the discretization error 
from the trapezoidal scheme has been significantly reduced, an additional error has 
appeared in the form of oscillations. Therefore, it was concluded that the trapezoidal 
scheme is not suitable for the solution of the governing equations resulting from these 

theories. 

The results shown in Fig. 5.8 are for the solution of the same test problem using a two- 
stage Gauss implicit Runge-Kutta finite difference scheme. This one-step scheme has an 
accuracy of order four and clearly provides a dramatic improvement for this test problem. 

In addition to superior accuracy as compared to the trapezoidal scheme, the Gauss scheme 
also offers improved stability for stiff ODEs due to its implicitness. This characteristic is 
desirable because DAEs are known to behave numerically similar to stiff differential 

equations. 51 
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5.7 Test problem results for w = 0.5 cm translation, solution by trapezoidal scheme 

solution with 320 intervals. 



Fig. 5.8 Test problem results for w = 0.5 cm translation, solution by two-stage Gauss 

scheme solution with 40 intervals. 
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This scheme is presented on pp. 210-216 of Ref. 52 and is summarized in Appendix 
B. In that form, this scheme is applicable to systems of first order linear ODEs. Therefore, 
the DAEs of the present theories were converted to ODEs by differentiating the algebraic 
equations, making them ODEs, and applying them in their original algebraic form as 
boundary (in our case initial) conditions. While this is unlikely to be an optimal method for 
the solution of these systems from the aspect of efficiency, it is sufficient for the 
verification and evaluation of the theories and is used in the remainder of this work for both 
the generalized plane deformation and axisymmetric models unless otherwise noted. 

5.4 Patch Tests 

Various test problems were solved by the structural models in order to verify their 
accuracy and to examine their applicability and limitations. The test problems examined are 
what will be referred to as patch tests because of their similarity to tests of that name used 
to verify the accuracy finite elements. The concept of a patch test is very simple: choose a 
problem for which the exact elasticity solution is known, divide it into layers of varying 
geometries, and compare the results obtained to the exact solution. Usually the exact 
solution is one with uniform strain and stress states. Using this approach, the effect of such 
factors as the layer distortion and thickness on the accuracy of the solution were examined. 

5.4.1 Generalized Plane Deformation Model 

Patch tests were used to evaluate the effect of the magnitude of rigid body 
displacements, steepness of the slope of interfaces between layers, the thickness of 
individual layers, and the number of layers in the model on the accuracy and reliability of 
solutions using the generalized plane deformation model. In each of these patch tests, the 
material was the same isotropic aluminum alloy used in the example of Sec. 5.3.2. In each 
case, the uniform thickness plate was subjected to a uniaxial strain of e a =0.1%, was 
constrained to £ w = 0, and was free to contract in the z-direction. This loading was chosen 
because it is characteristic of the type of loading applied to plates. The exact solution for 
this loading is <7„ = 76.281 MPa, = 23.647 MPa, and all other stresses zero. 

5.4. 1.1 Examination of the Effect of the Magnitude of Rigid Body 
Displacements 

The test problem used in examining the effect of the magnitude of the rigid body 
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displacements is shown in Fig. 5.9. The geometry is the same as that used in Sec. 5.3.2 
with only the boundary conditions different. The displacement in each coordinate direction 
is given the same amount of rigid body translation, A. The finite difference mesh in each 
case contained 20 equal sized intervals. 



Fig. 5.9 Patch test problem for rigid body displacements. 

The absolute value of the maximum error for each value of A was determined for the 
nonzero stresses o„ and o n as was the magnitude of the other stress components relative 
to the exact o„, since their exact values are zero. These results are presented in Table 5.2. 
It is clear that the finite difference solution by the two-stage Gauss implicit Runge-Kutta 
scheme does not suffer from difficulty in handling solutions with large amounts of rigid 
body translation as did the trapezoidal scheme. Even for the extreme case of A = 1000 cm , 
the errors all six stress components are, for all practical purposes, insignificant. 

Table 5.2 Errors resulting from rigid body translations. 


A (cm) 

<T„ (%) 

Oyy (% ) 






1.0 

0.0 

0.0 

0(1 0^) 

0(10“) 

O(10 7 ) 

0(10-“) 

10.0 

0.0 

0.0 

0(1 0^) 

0(1 0>°) 

0(10 7 ) 

O(10 10 ) 

100.0 

0.0013 

0.0 

0(1 0^) 

O(10- 9 ) 

0(10^) 

0(10 9 ) 

1000.0 

0.0315 

0.055 

1.6x10^ 

0(10*) 

0(10-5) 

0(10*) 
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5.4.1.2 Examination of the Effect of the Slope of the Interface Between 
Layers 

Due to the assumption that cosy? 1 0 for these models, a vertical layer interface 
(y = 90°) cannot be modeled. In this section the ability of the generalized plane deformation 
version to provide accurate solutions for very steep layer interfaces is determined. The test 
problem used for this is the same as that of Fig. 5.9 with one exception, the length of the 
plate will not be 10 cm. In order to increase the slope ( dzjdx ) of the internal interface, the 
length of the plate l is reduced, compressing the curvature of the internal interface into a 
shorter horizontal distance and therefore increasing the slope. The loading will remain at 
e =0.1% and A will be retained at 10 cm in order to examine the combined effect of rigid 
body displacements and a steep interfacial slope. Once again, the finite difference mesh will 
contain 20 intervals. 

The maximum errors in <7** and <7„, as well as the relative values of the other 
components were again calculated and are presented in Table 5.3. The maximum slope 
occurs at the mid-length of the plate and is also presented in the table. These results clearly 
show that this theory has the ability to model problems with very steep internal interfaces, 
likely beyond those required for the analysis of composite laminates. 


Table 5.3 Errors resulting from steep interfacial slopes. 


Max. 

\dzjdx\ 

On (%) 

Oyy (%) 

K/°«\ 

Kz/^| 

k*/*J 


1.125 

0.0 

0.0 

0( 1(H) 

0(1 0“ 10 ) 

0(1 0 7 ) 

0(1 0' 10 ) 

11.25 

0.0 

0.0 

0(10-*) 

0(io- 10 ) 

0(10*) 

0(1O 10 ) 

112.5 

0.0 

0.0 

0(10-*) 

0(1O- 10 ) 

0(1 0 7 ) 

0(1O 10 ) 

1125 

1.37 

1.43 

7.5X10- 4 

0(1 0-‘°) 

0(1(H) 

0(io- 8 ) 


5.4.1.3 Examination of the Effect of the Relative Thickness of Layers 


As mentioned in Sec. 3.3, the choice of displacement weighting functions in the two 
models presented in this work differ from those used by Pagano 44 - 45 in his original 
development of these types of theories. The reasons for the change centered on difficulties 
with numerical conditioning for problems with thin layers and for layers located at 
moderately large values of z relative to their thickness. Therefore, the first part of this 
section examines the degree to which the changes made were successful in overcoming 
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these difficulties. In the remainder of this section the accuracy of solutions using the higher 
order two-stage Gauss finite difference scheme for the same patch tests is examined. 

Two versions of the generalized plane deformation model will be compared in their 
reliability in solving problems with thin layers. The first was derived using the same 
approach as Pagano used for derivation of his model for involute bodies of revolution and 
will be referred to as Pagano’s approach. The second incorporates a series of modifications 
to the stress assumptions as well as the aforementioned change to the displacement 
weighting functions and will be referred to as the present approach. 

In his original development of the stress assumptions, Pagano did not scale the bubble 
shape functions by powers of the thickness. He also did not choose his cubic shape 
functions such that they had their inflection points within the layer. These changes to the 
stress shape functions were made in the present approach prior to the changes in the 
displacement weighting functions as an attempt to overcome the numerical conditioning 
problem encountered when solving the system, but were found to be insufficient. 

However, they did provide a side benefit in that they simplified the equilibrium equations 
and the expressions resulting from evaluation of the weighted compliance integral of Eq. 
(5-3). The changes made to the displacement weighting functions, on the other hand, were 
found to provide significant improvement in numerical conditioning as will be 
demonstrated. 

Both of these generalized plane deformation models were solved in their original DAE 
form using the trapezoidal one-step finite difference scheme. While this may not seem 
appropriate considering the inability of this solution approach to handle the test problem of 
Sec. 5.3.2, it is able to handle the present problem due to the absence of large 
displacements. For reasons that will become apparent later in this section, large 
displacements have little to no effect on the conditioning. Therefore, it was felt that the 
additional work involved in converting Pagano’s approach to a fully differential system to 
be solved by the higher order scheme is unnecessary for present purposes. 

The problem used to demonstrate the improvement in conditioning is shown in Fig. 
5.10. It differs from the previous patch test problems by the introduction of a third layer of 
thickness t whose mid-surface is defined by the same contour as the internal interface in the 
previous patch tests. Two geometric parameters are changed to examine their effect on the 
solutions by the two models. They are the thickness of the middle layer, which is varied by 
movement of the layer interfaces without moving the layer mid-surface, and the distance 
from the bottom of the plate to the x axis. This second parameter will be referred to as z 0 in 
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keeping with the notation used in presenting the model in Chapter 3 and is used to examine 
the effect of varying the magnitude of z within the layers relative to the thickness of the 
layers (Az). 



v = A v = A 

x z =0 i z =0 

Fig. 5.10 Patch test problem for examining relative thickness of a layer. 

The numerical conditioning of each case being referred to in this section is a measure of 
the reliability of the solutions obtained. The condition numbers that are presented are for the 
global finite difference matrix containing both the difference equation for each interval in 
the domain and the boundary conditions (see Eq. (5-15)). For the present class of 
problems, the boundary conditions are separated, leaving a banded matrix. The condition 
numbers presented here are the li condition numbers estimated by the IMSL library routine 
named DLFCRB 53 . This routine uses an algorithm developed by Cline et al. 54 for 
estimating the condition number of a matrix with minimal computational expense. 

The significance of the condition number in determining the reliability of the solution x 
to the general linear system 

Hx = b (5-21) 

can be expressed in terms of the sensitivity of the solution to perturbations in H and b. The 
condition number for H will be represented by k(H). In essence, the relative enror in x can 
be k(H) times the relative errors in H and b (see pp. 24—28 of Ref. 55). Therefore, a 
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matrix is considered to be ill-conditioned when the order of magnitude of the condition 
number approaches the number of significant digits used in the computations. For single 
precision fortran real numbers this is about seven and for double precision, which is used 
in the present study, it is about fifteen. Therefore, a condition number on the order of 10 15 
would indicate a poorly conditioned matrix and therefore an unreliable solution. It should 
be noted that while a poorly conditioned global finite difference matrix may produce reliable 
results for some right-hand-side vectors, this cannot be expected for all possible right- 
hand-side vectors as would occur with a well conditioned matrix. 

Having described how the numerical conditioning of a particular problem was 
examined, it becomes clear why the presence of large displacements does not affect the 
conditioning significantly. The only factor in the prescription of boundary conditions that 
affects the condition of the matrix H of Eq. (5-21) is the type of condition, i.e. traction or 
displacement. The magnitude of the prescribed quantities enters the equation through the 
vector b and therefore has no effect on the condition of H. 

The first set of results presented in Table 5.4 are for the case of A = 0 and Zo = 0 while 
the thickness of the center layer is varied. The solution is approximated by a uniform mesh 
of 41 grid points, giving 40 intervals. This choice of mesh refinement is somewhat 
arbitrary in this demonstration as it has been found to have little effect on the condition of 
the solution matrix. 


Table 5.4 Condition numbers for A = 0 and Zo = 0. 


t (cm) 

Pagano’s Formulation 

Present Formulation 

0.1 

3.5xl0 13 

3.7xl0 6 

0.01 

2.1xl0 20 

2.2xl0 8 

0.001 

7.1xl0 23 

1.4x10“ 

0.0001 

3.4x1 0 33 

2.0xl0 14 


From these results, it is clear that reducing the thickness of the center layer, at least 
relative to the other layers, causes a degradation in the condition of the finite difference 


solution matrix. It is also clear that the changes made to Pagano’s formulation in arriving at 
the present formulation significantly improve the conditioning. While Pagano’s formulation 
is on the verge of producing unreliable results for t = 0.1 cm, or 10% of the overall 
thickness of the plate, the present analysis is reliable past t = 0.001 cm or 0. 1% of the 
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overall thickness. 

The next set of results is for the case where z 0 = 100 cm, A = 0, and t is varied. For 
this case, the value of z within a layer has a very small range of values. Once again, a 
uniform mesh of 40 intervals was used in the approximation. The resulting condition 
numbers for this case with varying center layer thickness are presented in Table 5.5 for the 
two models. 


Table 5.5 Condition numbers for A = 0 and Zo = 100 cm. 


t (cm) 

Pagano’s Formulation 

Present Formulation 

0.1 

l.OxlO 17 

5.0xl0 6 

0.01 

5.3xl0 21 

2.6xl0 8 

0.001 

2.6x1 0 20 

1.4x10“ 

0.0001 

9.1xl0 29 

2.0xl0 14 


Compared to Table 5.4, these condition numbers reveal three noteworthy characteristics 
of these two models. First, the condition number for Pagano’s formulation with r = 0.1 cm 
deteriorates by three orders of magnitude to the point of producing unreliable solutions. 
Second, the present formulation of the model shows no degradation in the condition 
number relative to those shown in Table 5.4. Third, the trend of increasing condition 
number with decreasing t is reversed in going from t = 0.01 to t = 0.001 using Pagano s 
formulation. A re-examination of the two condition numbers in Table 5.4 associated with 
the same two thicknesses shows that the rate of increase in the condition numbers with 
thickness decrease is lower in that range for zo — 0 as well. This may be caused by a 
secondary influence on the condition numbers or possibly errors in the estimate of the 
condition numbers. However, because this anomaly occurs well into the range of condition 
numbers associated with ill-conditioned matrices, the matter was not examined further. 

The next set of results is for the case where t = 0. 1 cm, A = 0, and z 0 is varied. Once 
again, a uniform mesh of 40 intervals was used in the approximation. The resulting 
condition numbers for this case are presented in Table 5.6 for the two formulations of the 
model. 
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Table 5.6 Condition numbers for increasing distance from the x-axis, A = 0 and t = 0. 1 . 


z o(cm) 

Pagano’s Formulation 

Present Formulation 

0 

2.8xlO n 

3.7xl0 6 

1 

2.4xl0 14 

5.3xl0 6 

10 

8.3xl0 16 

4.4x1 0 6 

100 

l.OxlO 17 

5.0xl0 6 

1,000 

5.9xl0 18 

4.2xl0 6 

10,000 

5.3xl0 22 

4.4xl0 6 


These results reinforce the advantage of the present formulation over Pagano’s 
approach. It is clear that the present formulation’s condition number for the solution matrix 
is basically constant with respect to the z location of the domain. Pagano’s formulation, on 
the other hand, shows a degradation in conditioning with increasing z 0 . 

The remainder of this section will focus on the effect of the relative thickness of layers 
on the accuracy of the solutions obtained. The same patch tests used in examining the 
reliability of the method are solved using the higher order two-stage Gauss implicit Runge- 
Kutta finite difference scheme. As in Secs. 5.4. 1.1 and 5.4. 1.2, 20 finite difference 
intervals will be used and the maximum relative errors in and are presented. 

Table 5.7 presents the condition numbers and the accuracy of the stresses predicted by 
the model. Results for o yz and are omitted because the analysis gave values of zero for 
them, which is exact. Comparison of the condition numbers from the second column of 
Table 5.7 with those from the third column of Table 5.4 reveals a degradation in 
conditioning. This is due to a combination of two factors; the conversion of the system of 
DAEs to ODEs, and the use of the higher order Gauss finite difference scheme. A test 
application of the trapezoidal scheme to the fully differential system used with the Gauss 
scheme showed that the majority of the degradation was caused by the conversion to the 
fully differential system. The accuracy results show good agreement down to a thickness of 
the internal layer of t = 0. 00 1 cm , or 0. 1 % of the total thickness of the plate. 
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Table 5.7 Condition numbers and maximum errors for varying center layer thickness. 

r(cm) Condition# <r„ (%) (%) 

0.1 2.0xl0 10 0.0 0.0 0(10-®) 0( 1 0 7 ) 

0.01 2.9xlO n 0.0092 0.013 O(10 5 ) 0(1O- 5 ) 

0.001 2.1xl0 13 0.84 0.79 5.4x10^ 7.9x10^ 

0.000 1 2.7xl0 18 18 17 5.6xlQ- 3 9.0xlQ- 3 

An additional note should be made about the accuracy results presented in Table 5.7. 
Although the errors presented are large for the case of t = 0.0001 cm, these are maximum 
errors and are far from uniform over the domain. Fig. 5.11 is a plot of the errors in o xx 
along a series of parallel contours of constant x. From this plot, it is clear that the errors are 
mosdy confined to the thin central layer and that the solution is quite good everywhere else. 



x 

Fig. 5.11 Plot of errors in for t = 0.0001 cm case of Table 5.7. 


5.4.1.4 Examination of the Effect of Several Thin Layers 

The geometry of the modeled dropped-ply laminate of Fig. 5.1 demonstrates the need 
for this model to handle a group of very thin layers. Those layers used to model the 
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dropped plies cannot be discontinued and therefore are treated as a thin resin layer between 
the continuous sublaminates. This section includes such a feature in a patch test to examine 
its effect on the accuracy of the model. 

The problem used is an extension of the previous patch tests and is shown in Fig. 5. 12. 
The geometry differs from that of the previous section only in that it includes additional 
internal thin layers. Each of these internal layers has thickness 0.001 cm, or 0.1% of the 
thickness of the entire plate. As a group, these layers are centered within the plate, that is 
the layer geometry is again antisymmetric with respect to x = 5. The material and boundary 
conditions are the same as the previous section, i.e. aluminum subjected to 0.1% axial 
strain and constrained from contracting in the y-direction. 



Fig. 5.12 Patch test problem for examining several thin layers. 

Condition number and accuracy results for this patch test are presented in Table 5.8. 
These results were obtained using a finite difference mesh of 20 evenly spaced points. 
Once again, results for Oy, and are not presented because the theory correlated with the 
exact solution perfectly. The parameter N tabulated in the first column refers to the total 
number of layers in the model and therefore is two more than the number of interior layers. 
The case of A = 3 corresponds to the configuration in the previous section with 
t = 0.001 cm. 
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Table 5.8 Condition numbers and maximum errors for several thin layers of thickness 

t = 0.001 cm. 


N 

Condition # 

(%) 

o yy (%) 

ojdj 


3 

2.1xl0 13 

0.84 

0.79 

5.4x10^ 

7.9X10 -4 

4 

2.9xl0 14 

0.75 

6.0 

0.050 

8.6x10^ 

5 

6.7xl0 14 

0.25 

15 

0.014 

6.0X10 -4 

6 

l.lxlO 15 

0.27 

16 

0.051 

4.0x10^ 


The results shown in table 5.8 reveal a potential drawback to this model in the form of 
deteriorating condition numbers for increasing numbers of thin layers. In addition, the 


accuracy of the transverse and through the thickness stresses Oyy and o*, degrades with the 
number of thin layers. However, as noted previously, an ill conditioned matrix does not 
mean the solution will be poor, only more likely so. Therefore, care must be taken when 

evaluating solutions obtained with this model. 

As far as the accuracy of results for the stresses, all but and O a agree well with the 
exact solution. Figs. 5.13 and 5.14 are plots of these two stress components for the case 
N = 6 revealing once again that the errors in predicted stresses tend to be isolated within 
the extremely thin layers. While it would be desirable to obtain a solution that is uniformly 
good, knowing where stresses are likely to be in error allows values obtained there to be 
scrutinized more carefully. 
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5.4.2 Axisymmetric Model 

Patch tests were also used to examine the stability and accuracy of the axisymmetric 
model. The effects of the magnitude of rigid body displacements, the slope of interfaces 
between layers, and the thickness of individual layers in the model on the results obtained 
and their reliability were examined. Once again the material chosen was isotropic, with 
material properties of an aluminum alloy, E = 68.95 GPa and v = 0.31. The structure 
chosen for the tests was a uniform thickness cylinder of constant radius subjected to a 
uniform end extension of £„ = 0.1% and an internal pressure of p = 10 kPa. The cylinder 
was free to expand or contract radially along its entire length, including the ends. The exact 
elasticity solution for this loading is given by the expressions 56 


a„ - Ee r + vl 


y r N ~ r o 




(5-22) 


where r 0 and r N are the inside and outside radii respectively. 

5.4.2. 1 Examination of the Effect of the Magnitude of Rigid Body 
Displacements 

The test problem used in examining the effect of the magnitude of the rigid body 
displacements on solutions from the axisymmetric model is shown in Fig. 5.15. The 
geometry of the cross section is the same as that used in Sec. 5.4. 1 . 1 except that it now 
represents the cross section of a shell of revolution. The inside radius of the cylinder r 0 is 
200 cm, giving a radius to thickness ratio R/T of about 200, which was chosen because it 
represents a typical value for the unstiffened shells that are of interest. The amount of rigid 
body displacements in the x and 8 coordinate directions is denoted by A. For the 6 
direction, A is the circumferential displacement of the inner surface. There is no rigid body 
mode in the r coordinate direction because of the assumption of axisymmetry. The finite 
difference mesh in each case contained 20 equal sized intervals. 

The absolute value of the maximum error for each value of A was determined for the 
nonzero stresses and o n . For the other three stress components, their magnitude 

relative to the exact o a was determined, since their exact values are zero. These results are 
presented in Table 5.9. As in the generalized plane deformation model, the results are 
virtually unaffected by the presence of the rigid body displacements. The only components 
to show a degradation in accuracy as A was increased were and o x0 , but these errors 
were very small except for extremely large values of A. The errors appearing in the C7„ and 
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Gee stress components were due to the fact that only five significant digits of the stress 
variables were written to the output files from the analysis, which is why they do not 
change with increasing A. The error in o n is probably due to the approximate nature of this 
model and its numerical solution and is well within acceptable limits. 



x = Q Jt = 10 

u= A u = 0.01 + A 


v = (A/r 0 )r v = (A/r 0 )r 

T r =0 T r =0 


Fig. 5.15 Patch test problem for rigid body displacements in the axisymmetric model. 


Table 5.9 Errors resulting from rigid body displacements. 


A (cm) 

O a (%) 

Oee (%) 

O n (%) 


Ger/On 

\GjOn | 

\®x6 1® xx 


1.0 

7.8x10^ 

3.8xl0- 3 

9.0x1 0- 2 

4.9x1 0- 5 


9.5x10 - 

) 

10.0 

7.8x10^ 

3.8xl0- 3 

9.0xl0- 2 

4.9x1 0- 4 


9.5x10 4 

100.0 

7.8x10-* 

3.8xl0- 3 

9.0x10-2 

4.9X10- 3 


9.5x10-5 

1000.0 

7.8x10-* 

3.8xl0- 3 

9.0x10-2 

4.9x10-2 


9.5x10-2 


S.4.2.2 Examination of the Effect of the Slope of the Interface Between 
Layers 

In the same manner as was used to test the generalized plane deformation model, the 
effect of steepness in slope of layer interfaces was examined. The test problem used for this 
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is the same as that of Fig. 5.15 with A = 0 and with the length of the plate l successively 
reduced, thereby increasing the slope. The loading will remain the same as that applied in 
the previous patch test, end extension with e„ = 0.1% and internal pressure p = 10 kPa . 
Once again, the finite difference mesh will contain 20 intervals. 

The absolute value of the maximum errors in o a , OW, and o n were determined for each 
length l based on the exact elasticity solution and are presented in Table 5.10. The 
maximum slope occurs at the mid-length of the plate and is also presented in the table. The 
errors in the shear stresses were negligible and are therefore not presented. Clearly these 
results indicate that the axisymmetric model does not handle the presence of steep interfacial 
slopes very well for this loading case. It was found, however, that removal of the internal 
pressure loading leaving uniform axial extension virtually eliminated these errors. 


Table 5.10 Errors resulting from steep interfacial slopes. 


Max. | dr k /dx\ 

<x„ (%) 

Oee(%) 


0.1125 

7.5x10^ 

3.8xl0- 3 

9.0x10 2 

1.125 

2.1xl0- 3 

2.4xl0- 2 

25.1 

11.25 

6.5x10^ 

6.5xl0- 2 

80.9 

112.5 

6.4x1 0 3 

9.3X10 1 

263 

S.4.2.3 Examination of the Effect of the Relative Thickness of Layers 


Having established that the changes made to the displacement weighting functions for 
the generalized plane deformation model were successful in overcoming the difficulties in 
numerical conditioning associated with layers that are thin relative to their coordinate 
location, this section presents results for similar tests of the axisymmetric model. The 
problem used to examine the accuracy and reliability of the solutions is shown in Fig. 5. 16 
As in the generalized plane deformation model patch tests, a third layer of thickness t was 
introduced whose mid-surface is defined by the same contour (antisymmetric with respect 
to x = 5) as the internal interface in the previous patch tests. 
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Fig. 5.16 Patch test problem for examining relative thickness of a layer. 

The intent of this patch test was to verify that the axisymmetric model solved by the 
two-stage Gauss finite difference scheme produces results that are accurate and reliable for 
problems containing very thin layers in cylindrical shells with ratios of the radius to the 
thickness R/T of about 200. That is, geometries for which the thickness t is very small 
relative to the shell thickness T = 1.0 cm, and r 0 = 200 cm. However, in the first case 
examined, where t = 0.1 cm, unacceptable results were obtained. Presented in Table 5.1 1, 
these results show that while errors in and O m were less than 3%, the maximum errors 
in o„ exceeded 1000% and the condition number of the finite difference solution matrix 
exceeded 10 17 . As in the interface slope patch test, the three shear stress components were 
approximately zero, which is their exact value, and were therefore left out of this table. 

In order to establish the effect of reducing the thickness of the internal layer, a second 
case with the same R/T was examined with t reduced to 0.01 cm. The results are presented 
in the second row of Table 5.1 1 and, as expected based on the results from the generalized 
p lan e deformation model, they were worse than for the larger value of t. Both the condition 
number and the accuracy of the stresses deteriorated significantly. 

To determine the effect of varying the shell R/T ratio, the value of r 0 was reduced to 20 
cm and t was returned to 0.1 cm while all other dimensions were left unchanged, giving an 
R/T of about 20. The results for this case are given on the third row of Table 5.11. 
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Compared to the first row of the table, the solution matrix remained poorly conditioned and 
the accuracy of the solution improved somewhat. While the errors in O a and o ee were 
acceptable at less than 1%, the error in was again much larger at 27%. 


Table 5.11 Condition numbers and maximum errors for varying center layer thickness 

and radius. 


t (cm) 

r 0 (cm) 

Condition # 

<*„(%) 

Oee(%) 

<x„(%) 

0.1 

200.0 

2.0xl0 17 

0.093 

2.59 

1480 

0.01 

200.0 

1.3X10 20 

232 

2560 

2.20x1 0 5 

0.1 

20.0 

1.7xl0 16 

0(10-*) 

0.146 

27.2 

1.0 

200.0 

3.2xl0 14 

0(1O- 5 ) 

3.11 

78.4 


Finally, all of the radial dimensions were scaled by a factor of ten from the previous 


case, leaving a cylinder of inside radius r 0 — 200 cm , total thickness T — 1 0 cm with a 1 cm 
thick inner layer. These results are given in the last line of Table 5.1 1 and actually show a 
deterioration in accuracy despite an improvement in the condition number. This increase in 
the errors may be due to the increase in interfacial slopes that this scaling produces, which 
was shown previously to affect the accuracy in the presence of an internal pressure loading. 
Removal of the internal pressure, as before, virtually eliminated these errors. 

5.4.2.4 Comparison of Condition Numbers to Pagano’s Original 
Axisymmetric Formulation 

Although this reformulation of Pagano’s original axisymmetric model still suffers from 
poorly conditioned finite difference solution matrices for the types of geometry necessary to 
analyze laminated cylindrical shells with dropped plies, it was successful in significandy 
improving the situation. This section presents a comparison of finite difference solution 
matrix condition numbers between the present axisymmetric formulation and Pagano’s 
original model. The reader is referred to Sec. 5.4. 1.3 for a brief explanation of the changes 
made to Pagano’s approach or to Chapter 4 for a more thorough treatment 

As in Sec. 5.4. 1.3, the comparison of conditions numbers was done using the one-step 
trapezoidal finite difference scheme as applied to the original DAE form of the governing 
equations. The problem used to demonstrate the effect of the reformulation on the 
numerical solution conditioning has the same two-layer geometry as the patch test problem 
used for examination of the effect of rigid body displacements and is shown in Fig. 5.15. 
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However, the total cylinder thickness is not 1 cm for all of these cases. The finite difference 
mesh contained 80 intervals in each case. 

The condition number results for this comparison are presented in Table 5.12. As the 
accuracy of the results were bad for all cases due to the use of the trapezoidal scheme, they 
are not presented. The first and second cases differ only in that the radius of the cylinder is 
reduced from 200 to 20 cm. The third case is a radially scaled version of the second by a 
factor of 10. In a similar fashion to the results presented in Sec. 5.4. 1.3 for the generalized 
plane deformation model, the condition numbers for the present version of the model are 
much lower than those for Pagano’s formulation. They also tend to be influenced much 
less by the thickness of the cylinder and its radius. Another similarity between the condition 
number results from the axi symmetric model and the generalized plane deformation model 
is that the change in solution method from the trapezoidal scheme applied to the DAEs to 
the more accurate Gauss finite difference scheme applied to a fully differential system 
causes a significant degradation in the condition number. For example, the condition 
numbers for the rigid body displacement patch test were on the order of 10 14 while the 
same geometry solved by the trapezoidal scheme, the first row of Table 5.12, had a 
condition number in the 10 9 range. 


Table 5.12 Condition numbers for Pagano’s original formulation and the present 

formulation. 


T( cm) 

r 0 (cm) 

Pagano’s Formulation 

Present Formulation 

1.0 

200.0 

7.6x1 0 26 

3.6x1 0 9 

1.0 

20.0 

3.6xl0 29 

2.2xl0 9 

10.0 

200.0 

1.4x1 0 23 

1.7xl0 9 


S.4.2.5 Comments Regarding the Patch Test results for the Axisymmetric 
Model 

As the results for the case where the thinnest layer was 10% of the total thickness with 
an RJT of 20 or 200 were unsatisfactory, modeling of resin layers on the order of 10% the 
thickness of a single ply in a composite laminate with a similar R!T does not appear to be 
feasible for an internal pressure loading. Despite the improvement in the conditioning of the 
numerical approximation attained through the modifications found to be successful in the 
generalized plane deformation, solution stability and accuracy are still deficient in the 
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present formulation and its solution by finite differences. The cause for this difference in 
the level of success likely lies in the inherent differences between plates and shells. 

Variation of r 0 , while somewhat similar to the variation of Zo for the flat plate of the 
generalized plane deformation model, affects more than the numerical representation of the 
problem and its solution, it also affects the solution itself. That is, varying r 0 while keeping 
everything else the same changes the curvature of the shell and therefore, in general, the 
stress and displacement fields. This effect is magnified by the presence of an internal 
pressure loading. On the other hand, varying zo only changed the way the generalized plane 
deformation model represented the solution to what was essentially the same problem. The 
modifications to Pagano’s approach simply minimized the degree to which varying z<, 
changed that representation. 

The main purpose of developing the axisymmetric model was for the analysis of 
laminated shells of revolution with dropped plies subjected to internal pressure. While the 
present formulation produces much better patch test results without internal pressure than 
with, the errors occurring with internal pressure are disconcerting. Therefore it was decided 
further application of this model, even to problems without internal pressure, should not be 
trusted without a better understanding of the cause of these errors. 

5.5 Comparison with Other Solutions 

While the patch tests of the previous section are useful for establishing some confidence 
in the generalized plane deformation model, they do not go far enough. Because their 
premise requires an exact solution to compare their results to, they are limited to very 
simple geometries. Problems involving multiple constituents for which exact elasticity 
solutions are available do not exist. Therefore, in order to verify the ability of these models 
to predict the stresses in composite laminates, comparisons must be made with established 
approximate solutions. This section presents results for the tensile coupon free edge 
problem and also for the problem of a dropped-ply laminate under compression, both of 
which are solved using the generalized plane deformation model. 

5.5.1 Tensile Coupon Free Edge Problem 

The first verification of the structural model as applied to laminates was a comparison to 
the benchmark problem of the interlaminar stress response near the straight free edge of a 
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tensile coupon. This was chosen as a baseline case because of the simple geometry and the 
availability of solutions. A four ply, symmetric laminate with [0/90] s layup is the test case 
selected (A 0°-ply is parallel to the y-axis). Each ply has uniform thickness h, the width of 
the laminate is 2b, and b = 8 h. See Fig. 5.17. The load is a prescribed normal strain e in 
the y-direction, so C = e in Eqs. (3-1) and (3-2). All four lateral surfaces are traction free. 
We compare our results to those published by Pagano 44 for the interlaminar stress 
distributions along the x-axis in the 0/90 interface (z = h). Ply material properties are 

E n = 138 GPa (20 Msi), E 22 = = 14.5 GPa (2.1 Msi) 

G l2 = G 13 = G 23 = 5. 86 GPa (0. 85 Msi) 
v 12 = V !3 = v 23 =0.21 

Because of problem symmetry about the x-y and y-z planes, the solution domain is reduced 
to the quadrant 0 <x<b and 0 < z < 2/i in the y = 0 plane. Each ply is divided into three 
mathematical layers, so N = 6. 


A 



Fig. 5.17 Four ply laminate subject to uniform extension on ends y = constant. 

The interlaminar normal stress distribution is shown in Fig. 5.18 and the distribution of 
the interlaminar shear stress, z xz , is shown in Fig. 5.19. It is clear from these plots that 
results of the present model compare very well to those presented by Pagano. A good 
comparison was expected for this problem since the stress assumptions in our structural 
model for straight and uniform thickness layers essentially reduce to those utilized by 
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o z (x,h)/10 6 e, kPa 


Pagano in Ref. 44. It is worth noting that Pagano compared his results to those of Wang 
and Crossman 57 and found excellent agreement. 



Fig. 5.18 Distribution of the interlaminar normal stress <7 Z in the 0/90 interface. 
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Fig. 5.19 Distribution of the interlaminar shear stress x xz in the 0/90 interface. 

5.5.2 Dropped-Ply Laminate Under Uniaxial Compression 

The second test case selected for the structural model as applied to laminates is a 
comparison to the results of Curry et al. 21 for the interlaminar stress response of a 
dropped-ply laminate subjected to axial compression. This problem was solved by a 
displacement-based finite element analysis in Ref. 21, and the finite element analysis was a 
part of an effort to correlate the prediction of delamination initiation with test results. Thus, 
the problem definition (laminate configuration, boundary conditions, etc.) were dictated by 
test conditions. 

The laminate analyzed had a 20-ply thick section with the layup [(±45/0/90) s /0 2 ] s from 
which the center four 0° plies were dropped leaving a 16-ply quasi-isotropic [±45/0/90]^ 
layup in the thin section. The thick portion is 2.79 mm (0.1 1 in.) thick and the thin portion 
is 2.24 mm (0.088 in.) thick. Each ply was modeled with one element through the 
thickness with the exception of the two adjacent 90° plies which were modeled together by 
a single element through the thickness. The portion of the mesh used to model the transition 
region is shown in Fig. 5.20. 
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Fig. 5.20 Finite element mesh for the dropped-ply laminate modeled in Ref. 21. 

This laminate was modeled with six mathematical layers as shown in Fig. 5.21. The 
two middle layers are given the material properties of 0° plies up to the drop after which 
they arc treated as resin. While the two outer sublaminates appear to rejoin in Fig. 5.21, 
they are actually separated by the two thin resin layers, each 6.35 |im (0.00025 in.) thick, 
still present in the model. Note that the layers closer to the stress concentration are thinner 
than the top and bottom layers. The layups for the four continuous layers are, from bottom 
to top, [±45/0/90 2 /0], [+45], [±45] and [0/90 2 /0/+45]. These layers, being composed of 
multiple plies of different orientations, are replaced by a single layer with effective 
compliances found using the homogenization procedure described in Sec. 5.2.2 and 
Appendix A. 

The finite element model contained 2,905 nodes with three degrees of freedom per node 
for a total of 8,715 degrees of freedom. The six layer model was solved with 58 finite 
difference mesh points for a total of 8,004 degrees of freedom. No comparison is made in 
computation time because the two models were solved on different machmes. However, 
because of the high cost of computing the difference equations using the two-stage Gauss 
scheme and the use of a banded matrix solver rather than a more efficient skyUne solver for 
the six layer model, it is estimated that the six layer model would require at least twice as 

much time. 
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Fig. 5.21 Ply drop-off as modeled using present analysis. 

The material properties for the individual plies are those of AS4/3502 graphite/epoxy: 

E n = 128 GPa (18.5 Msi), = E 33 = 11.3GPa (1.64 Msi) 

G n =G l3 =6.0 GPa (0.87 Msi), G 23 = 3.38 GPa (0.49 Msi) 
v 12 = v 13 = 0.3, v 23 =0.35 

and for the neat resin: 


£ = 3.45 GPa (0.5 Msi) 
v = 0.41 

The loading for the case presented is uniaxial compression and the results are normalized 
by the average normal stress in the x-direction in the thin section of the laminate, denoted 
by o app . The boundary conditions at the ends are displacement prescribed conditions 
derived from a plate analysis of the test specimen. 

The results of the analysis are shown in Figs. 5.22 through 5.25 and are taken along 
contours passing through the Gauss points of the finite element analysis. The terms lower 
and upper interface refer to the interfaces between the four 0° plies and the quasi-isotropic 
sublaminates. The stresses for the lower interface are taken along the contour passing 
through the Gauss points at the top of the elements modeling the 45° lamina just below the 
four 0° plies. The stresses for the upper interface are taken along the contour passing 
through the Gauss points at the bottom of the elements modeling the 45° lamina just above 
the four 0° plies. The stresses are taken at the Gauss points in order to obtain the most 
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accurate values from the fmite element analysis. The t-n notation refers to contour 
following coordinates with n normal to the contour and 7 tangent to it. The positive senses 
for them are shown in Figs. 5.20 and 5.21. The longitudinal location in these plots is 
normalized by the thickness of the four dropped plies, denoted by to (to = 0.5588 mm in 
this example). 
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Fig. 5.25 Comparison to finite element model of Ref. 21 for o m along upper interface. 


It is clear from Figs. 5.22-5.25 that the stresses predicted by the two approaches are in 
good agreement. While the two analyses do not agree as well as those in the previous 
section, this is understandable due to the coarseness of the finite element mesh used by 
Curry et al. 21 
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Chapter 6: Parametric Studies 


Having established the accuracy of the generalized plane deformation model, this 
chapter presents the results of parametric studies of dropped-ply laminates. A series of 
parametric studies were performed in order to examine the effects that various factors have 
on the stress concentration, and hence the potential for a reduction in strength, of laminates 
with dropped internal plies. The major factors that are felt to contribute to the stress 
concentration in an asymmetric dropped-ply laminate are the eccentricity of the middle 
surface through the ply drop-off and the magnitude of the stiffness discontinuity. 

The effect of eccentricity is examined first as its effect on the stresses can easily be 
isolated by using the same layups for varying amounts of eccentricity. The magnitude of 
the stiffness discontinuity for a laminate composed of one fiber/matrix combination, on the 
other hand, cannot be changed without also changing the orientation of plies or the relative 
thickness between the continuous and dropped sublaminates. Rather than altering the 
geometry of the transition region by changing the number of dropped plies, the approach of 
altering the orientation of the dropped plies was chosen. Finally, the effect of the inclusion 
of a soft insert ahead of stiff dropped plies was evaluated. 

6.1 Effect of Eccentricity of the Middle Surface 

When a laminate contains an asymmetric ply drop-off, the middle surface is not a 
straight line (in the plane of the cross section) but instead contains an offset at the drop-off. 
This eccentricity in the middle surface would appear to be a factor in the strength 
degradation that a ply drop-off can cause in a laminate due to the bending or twisting 
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moment that may be induced. This study was done to determine the effect that this 
eccentricity has on the interlaminar stresses present at a ply drop-off. 

Two geometries of ply drop-offs are compared, one that is asymmetric and therefore 
has an eccentric middle surface, and the other that is symmetric with no eccentricity (see 
Fig. 6.1). These two geometries were examined with three different layups for the dropped 
plies and under three different loadings for a total of eighteen different cases. 

t z 



t ^ 


[±45/0/90] s 

< 1 ' 

> 

j> 4 dropped plies | — — *i 


[±45/0/90] s 



X 


(b) Without eccentricity. 


Fig. 6.1 Laminates with and without eccentricity in the middle surface. 

In each case, the two continuous sublaminates have the same eight ply quasi-isotropic 
[±45 / 0 / 90] s layup. The dropped sublaminate contains four plies with the layups [0 4 ] T , 
[90J X , or [±45] s . The material properties used are those for AS4/3502 graphite/epoxy and 
are given in Sec. 5.5.2. Each ply has thickness 0.1395 mm (0.0055 in.) giving a thickness 
of 0.558 mm (0.022 in.) for the dropped plies, 2.794 mm (0.1 10 in.) for the thick portion, 
and 2.248 mm (0.0885 in.) for the thin, which also contains a 0.0127 mm (0.0005 in.) 
thick resin layer between the continuous sublaminates. The eccentricity is, therefore, either 
zero or 0.273 mm (two ply thicknesses minus half the resin layer thickness). The length of 
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the thick and thin portions is 25.4 mm (1.0 in.) for each and the transition section is 2.24 
mm (0.088 in.) long for a total laminate length of 53.04 mm (2.088 in.). The origin of the 
coordinate system for the study is located in the plane of the bottom of the thick portion of 
the laminate at the end of the dropped plies, i.e. the beginning of the thickness taper (see 
Fig. 6.1). 

The model used in this analysis consists of six mathematical layers with the same 
distribution as those used in Sec. 5.5.2. That is, in moving from one surface to the other of 
the thick section, the layers are [+45 / 0 / 90 2 / 0] , [+45], [±0], [+0], and 
[0 / 90 2 10/ +45] . From the drop into the thin section, the central two layers are modeled 
as resin except for the 90° case. Experience in manufacturing specimens with 90° dropped 
plies has shown that these fibers tend to migrate into the area occupied by resin in the [0 4 ] T 
and [±45] s layups. For this case, the triangular shaped resin region is modeled as 90° 
material with the origin remaining at the thick end of the taper. 

The three loading cases are longitudinal compression (along the x-axis), transverse 
tension (along the y-axis), and shear in the x-y plane. The surfaces of the laminates are 
considered to be traction-free in all three cases. The boundary conditions at the left (thick) 
end (jc = -25. 4 mm ) are also the same in all three cases, with u = 0, v = 0 , and t z = 0 . 
These conditions were chosen to avoid a boundary layer at the end by simulating conditions 
at a far field location. Each laminate is constrained from rigid body translation in the z- 
direction through the prescription of w = 0 at the lower surface of a very short (0.254 mm) 
segment at the left end. The other rigid body modes are excluded through the displacement 
boundary conditions. The longitudinal compression and in-plane shear loads are applied 
through prescription of N x and respectively at the right end (jc = 27.64 mm), while 
restraining the u displacement to a uniform value through the thickness. The transverse 
tension load is applied through prescription of the uniform strain £y out of plane 
deformation allowed for by the generalized plane deformation formulation, i.e. C of Eqs. 
(3-1) and (3-2) is prescribed. The boundary conditions and out of plane deformations used 
for all three loading cases are presented fully in Table 6. 1 . It is worth noting that the 
conditions applied on the upper and lower surfaces of the laminates are not boundary 
conditions in that they are actually applied over the entire one dimensional mathematical 
domain and are therefore part of the field equations. 
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Table 6.1 Boundary conditions and out of plane deformations for eccentricity study. 


Load Case 

BCs at x = -25.4 mm 

BCs at x = 27.64 mm 

Out of Plane Defs. 


u = 0 0 

N = - 1 kN / m 

A = 0.0 

longitudinal 

v = 0.0 

X 

u = uniform 

O 

O 

II 

CQ 

compression 

T = 0.0 

< 

II 

O 

o 

O 

O 

II 

O 


z 

T z = 0.0 

0 = 0.0 


u = 0.0 

o 

O 

II 

o 

o 

II 

transverse 

v = 0.0 

u = uniform 

to 

II 

p 

o 

tension 

T z =0.0 

v = 0.0 

C = 0.001 



o 

o 

II 

0 = 0.0 

in-plane shear 

u = 0.0 
v = 0.0 

N,= 0.0 

u = uniform 

A = 0.0 
0 = 0.0 


T z =0.0 

N xy = 1 kN / m 
v = uniform 

C = 0.0 
0 = 0.0 



o 

o 

II 



A brief note is in order to explain the x-direction conditions at x = 27.64 mm for all 


three load cases as well as the y-direction condition for the shear loading case. The 
combination of a uniform u or v displacement and N x or N„ prescription was achieved as 
follows. The x-direction boundary conditions of Eq. (3-52) require 2 N equations for the 
o[ k \ U a \ and U (k) unknowns. Constraining the u displacement to a uniform value 
through the thickness, U , gives the expressions for the weighted displacements 

V {k) = \/2U k = l ^ N (6-1) 

U (k) =l/2U 

from Eq. (3-14). Elimination of 0 from these 2N equations leaves the 2N-1 equations 

U (k) = U (k) k = l^>N 

£/<*> = k = l-+N-l 

The final equation is obtained through the definition for N x , 

N x = f N a„dz ( 6 ‘ 3 ) 

J Z 0 

which, when integrated through the thickness, leaves 
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(6-4) 


V 2 X( CT " , + <7 1 < * , )( Z * -Zk-l) = Nx 

k=\ 

due to the linearity of G a within layers. The tilde, as usual, indicates the prescribed value 
of N x . The same procedure was used in the y-direction for V ( *\ V tk) , & 6 \\ and cr^' i n 
terms of N . 

Since the primary focus of this study is the effect of eccentricity on the tendency of 
laminates with dropped plies to delaminate, a measure of this tendency is necessary. Rather 
than looking at each interlaminar stress component separately, it is desirable to have an 
index that accounts for the combined effect of the different components. An index similar to 
the Quadratic Delamination Criterion proposed by Brewer and Lagace 29 was chosen. This 
delamination index is defined as 



where Z s \ 2 s1 , and Z T are the allowable interlaminar stresses. For AS4/3502, they are 

Z sx = Z S2 = 93.08 MPa (13.5 ksi) 

Z T = 51.99 MPa (7.54 ksi) 

The onset of delamination is likely to occur when F > 1. However, F = 0.5 does not 
correspond to half of the delamination load for a linear analysis due to the quadratic nature 
of the index. Therefore, it was decided to examine the distribution of VF, which will be 
referred to as the delamination fraction. For a proportional loading, a base load resulting in 
a maximum value of VF can be increased by the factor 1 /Vf before delamination initiates. 
Hence, the larger the value of VF , the more likely delamination is to occur. 

The delamination fraction was examined along the two critical interfaces, i.e. the top of 
the lower continuous sublaminate and the bottom of the upper continuous sublaminate. 
These two interfaces will be referred to as the lower and upper interfaces respectively. For 
the non-eccentric cases, the distribution of VF is the same for the two interfaces due to 
symmetry. 

Results for the case of four 0° plies dropped from the laminate and subjected to 
longitudinal compression are presented here. Figs. 6.2 and 6.3 reveal that eccentricity has 
essentially no effect on the peak value of VF along either interface for this case (t d is once 
again the thickness of the dropped plies). There is a significant difference in the 
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interlaminar stresses in the transition region (0 < x / t d < 4), however these stresses are 
well below the peak values at the end of the dropped plies and are therefore not considered 

critical. 



Fig. 6.2 Delamination fraction along the lower interface for four 0” plies dropped loaded 

in longitudinal compression. 


It is informative to examine the interlaminar stress components separately to determine 
what each one is contributing to the delamination fraction. Figs. 6.4 and 6.5 show the 
distribution of o nn and a m (cr yn is approximately zero) along the upper interface (This 
loading represent an average compressive o„ stress in the thin section of the laminate of 
445 kPa.) The maximum value of the delamination fraction coincides with the maximum 
values of both the interlaminar normal and shear stresses, which also are not significantly 
influenced by eccentricity. This indicates that for this layup and loading, the stiffness 
discontinuity influences the interlaminar stresses much more than the eccentricity. As for 
the relative contribution of the two stress components, the peak interlaminar normal tensile 
value of 97 kPa is 0.18% of Z r while the peak interlaminar shear is -150 kPa, or 0.16% of 
Z s/ . This suggests that these components will have a comparable influence on the 
delamination of these laminates. 
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(kPa) 



Fig. 6.3 Delamination fraction along the upper interface for four 0° plies dropped loaded 

in longitudinal compression. 



Fig. 6.4 <j m along the upper interface for four 0° plies dropped loaded in longitudinal 

compression. 
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50 



x/t d 


Fig. 6.5 a m along the upper interface for four 0° plies dropped loaded in longitudinal 

compression. 

Up until about x = 2t d , the delamination fraction along the upper interface is very 
similar for the eccentric and non-eccentric cases. Each has a maximum at x = 0 with a 
second peak at about x = t d (the plots of the interlaminar stresses show that this second 
peak is caused by a second peak in the tensile interlaminar normal stress). The plots of the 
interlaminar stresses also closely agree up until about x = t d . For t d < x < 2 t d , the normal 
stress is less with eccentricity than without while the shear has a larger magnitude. These 
opposite trends tend to cancel leaving the failure fraction for the two cases in close 
agreement 

The delamination fraction is affected by the addition of eccentricity in the region 
2 t d <x< 4 t d where an additional peak arises. Examination of Figs. 6.4 and 6.5 clearly 
reveals that this second peak is due to an increase in the interlaminar shear stress in that 
region. A plot of the bending moment resultant in Fig. 6.6 shows a reversal in the bending 
moment through the transition region for the case with eccentricity (there is no bending 
moment for the case without eccentricity). This reversal in bending moment introduces the 
additional negative shear that produces the additional peak in the delamination fraction. 
While this additional peak is an interesting consequence of the presence of eccentricity, its 
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magnitude is still slightly less than the peak in the transition region that occurs without 
eccentricity and much less than the maximum that is located at the end of the dropped plies. 



Fig. 6.6 M x distribution for four 0° plies dropped loaded in longitudinal compression, 

with eccentricity. 

The results for all of the load cases and layups are presented in Table 6.2. The values of 
VF for load case b, transverse tension applied through prescription of £y, are scaled such 
that they correspond to N — 1 kN / m in the thin section. This was done in order to make 
the values of y[¥ correspond to comparable loads between the different load cases. 

For the longitudinal compression loading, the [90 4 ] T dropped layup experiences the 
largest change in the delamination fraction due to eccentricity, a 15% increase from 
1.22xl(H to 1.4xlCH. However, experimental results presented by Curry et al. 21 indicate 
that this configuration does not tend to fail by delamination under longitudinal 
compression. This layup having the smallest change in longitudinal stiffness and also the 
smallest values of Vf again indicates the larger role that the stiffness discontinuity plays in 
the delamination of dropped-ply laminates. The third layup ([±45] s dropped) loaded in 
longitudinal compression only experiences an increase of 3% in VF . 
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Table 6.2 Results for eccentricity study ( t lhx „ is the thickness of the thin section). 


Load Case 

Layup 

Eccentricity 

Peak VF 

x-location 

Interface 

Nx = -1 kN/m 

[0 4 ] T 

0.0 

2.25x10 3 

-0.077r d 

both 


0.121 t lhin 

2.25xl0- 3 

-0.077 

both 


[90 4 ] t 

0.0 

1.22xl0 3 

o 

oo 

o 

both 


0. 121 tlhin 

1.40x1 0- 3 

0.60t 4 

upper 


[±45] s 

0.0 

1.63xl0- 3 

-0.071 t d 

both 


0.121r lAj „ 

1.68xl0- 3 

0.0 

lower 

£y = 0.001 

[0 4 ] T 

0.0 

0.54X10- 3 

1 

O 

both 

Q.\2lt,hi n 

0.56xl0- 3 

-0.l5t d 

upper 


[90 4 ] t 

0.0 

0.09X10 3 

3.20t d 

both 


0*121 tfhin 

O.lOxlO- 3 

3.20t d 

upper 


[±45] s 

0.0 

0.42xl0- 3 

0.0 

both 


0.\2\t lhl „ 

0.45xl0- 3 

0.0 

lower 

Nxy = 1 kN/m 

[0J T 

0.0 

1.14xl0- 3 

0.0 

both 

0. 121 tfhin 

1.21xl0- 3 

0.0 

lower 


[90 4 ] t 

0.0 

0.31xl0- 3 

2A0t d 

both 


0. 121 t th i„ 

0.51xl0- 3 

2.30t d 

“PP er 


[±45] s 

0.0 

3.54xl0- 3 

0.0 

both 


0.12 \tfhin 

3.65xl0- 3 

0.0 

lower 


The values of 4f for the transverse tension loading are all well below the values for 
the [90 4 ] t dropped layup loaded in longitudinal compression, which again are below that 
required for a delamination failure at this level of load. Therefore, none of these three 


layups subjected to transverse tension would be expected to delaminate. Even so, 
eccentricity at most produces an increase in VF of 1 1% for layup 2. 

The shear loaded case also has its largest increase in Vf due to eccentricity occurring 
for the [90 4 ] t drop. This is also the case with the smallest value for the delamination 
fraction, too small to induce a delamination failure. The [±45] s drop layup, the stiffest with 
respect to shear, has the largest values of Vf for any of the cases, 3.65x1 0~ 3 for the 
eccentric laminate. However, this represents an increase of only 3% over the non-eccentric 

case. 
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Therefore, it appears that eccentricity does not significantly affect these laminates ability 
to resist delamination initiation when loaded by longitudinal compression, transverse 
tension, or in-plane shear. This does not imply that these laminates will not delaminate. 
Only that the presence of eccentricity of the middle surface does not appear to make this 
mode of failure more likely than for laminates without eccentricity. 

6.2 Effect of the Magnitude of the Stiffness 
Discontinuity 

Clearly, the results of the study examining the effect of the eccentricity of the ply drop- 
off indicate that the stiffness discontinuity has a much greater influence on the magnitude of 
the interlaminar stresses than the eccentricity of the dropoff. The purpose of the present 
study is to examine the relationship between the stiffness of the dropped plies and the 
resulting interlaminar stresses. 

The geometry of the laminate examined is the same as the eccentric dropped-ply 
laminate of the eccentricity study and is shown in Fig. 6.7. Four plies are dropped from 
between two eight ply quasi-isotropic [±45 / 0 / 90] s sublaminates with the bottom surface 
of the laminate flat. The stiffness of the four dropped plies is varied by altering the angle 6 
in the balanced angle-ply layup [±0] s . Five orientations for the angle 0 were examined: 0°, 
30°, 45°, 60°, and 90°. These layups were subjected to longitudinal compression and in- 
plane shear loadings for a total of ten different cases. The transverse tension loading 
examined in the eccentricity study was not repeated here because of the relatively low 
interlaminar stresses found for that loading. 

t ^ 



Fig. 6.7 Laminate used to examine effect of the stiffness of the dropped plies on the 

interlaminar stresses. 
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The material properties used are again those of AS4/3502 graphite/epoxy given in Sec. 
5.5.2. The dimensions of the laminate and the discretization used in modeling it are the 
same as those of the eccentricity study. The modeling of the 6 = 90 case, however, does 
represent a change worth noting. In the eccentricity study, the 90 fibers were assumed to 
migrate into the region occupied by resin for the other orientations. The present study , 
being more interested in determining the effect of the magnitude of the stiffness 
discontinuity, uses the same location for the material discontinuity (jc = 0) for 6 = 90 as 
the other layups for the sake of consistency and retains the triangular shaped resin region. 
The surface conditions are again traction free and the boundary conditions are the same as 
those given in Table 6.1 for the longitudinal compression and in-plane shear loadings of the 
eccentricity study. 

The measure of the severity of the interlaminar stresses used in evaluating these layups 
remains the delamination fraction (VF ). Figs. 6.8 through 6.11 are plots of the 
distribution of VF for the two loading cases. As in the eccentricity study, the lower and 
upper interfaces referred to are the interfaces between the dropped plies and the lower and 
upper continuous sublaminates, respectively. 

Figs. 6.8 and 6.9 show the distribution of the delamination fraction for the five 
different dropped-ply layups along the lower and upper interfaces for the longitudinal 
compression load case. Clearly, the peak value of the delamination fraction increases with 
increasing dropped sublaminate longitudinal stiffness, as expected. In addition, the decay 
length into the thick section required for the induced interlaminar stresses to vanish also 
increases. The interlaminar stresses within the transition region (0 < xjt d < 4) do not vary 
significantly with the stiffness of the dropped plies. 

Similar results are obtained for the in-plane shear loaded case shown in Figs. 6.10 and 
6.1 1 with the in-plane shear stiffness of the dropped plies as the major factor influencing 
the magnitude of VF. That is, the 0 = 45° case has the largest peak in VF as well as the 
greatest in-plane shear stiffness. As the value of 6 increases or decreases from 45 , the in- 
plane stiffness of the dropped plies decreases as does the peak v/F. An interesting 
difference between the response under in-plane shear and longitudinal compression 
loadings is that the decay length the interlaminar stresses is shorter for the shear loading. 
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-10 -8 -6 -4 -2 0 2 4 6 8 10 

x/t d 


Fig. 6.8 Delamination fraction along the lower interface for different values of 6 of the 
dropped sublaminate [±0] s (longitudinal compression loading). 



Fig. 6.9 Delamination fraction along the upper interface for different values of 8 of the 
dropped sublaminate [±0] s (longitudinal compression loading). 


Parametric Studies 


121 




gME 


Fig. 6.10 Delamination fraction along the lower interface for different values of 6 o e 
dropped sublaminate [i 0] s (in-plane shear loading). 


0 

= 0 ° 

0 

= 30 ' 

0 

= 45 ' 

0 

= 60 

0 

= 90 


Fig. 6.11 Delamination fraction along the upper interface for different values of d of the 
dropped sub laminate [i 0] s (in-plane shear loading). 
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The numerical values for the peak delamination fraction in each of the cases as well as 
the x-location and the interface where the peak value occurs is presented in Table 6.3. 


Table 6.3 Results for study of the effect of stiffness of the dropped plies ±6 (r, hm is the 

thickness of the thin section). 


Load Case 

e 

Peak 4f 

x-location 

Interface 

Nx = -1 kN/m 

0° 

2.25xl0- 3 

-0.077^ 

both 

o 

o 

1.96x10-’ 

1 

upper 

45° 

1.68x10’ 

0.0 

lower 

60° 

1.52x10’ 

0.0 

lower 

90° 

1.43x10’ 

0.0 

lower 

Nxy = 1 kN/m 

0° 

1.21x10’ 

0.0 

lower 

U> 

o 

o 

3.14x10’ 

0.0 

lower 

45° 

3.64x10-’ 

0.0 

lower 

60° 

3.44x10’ 

0.0 

lower 

90° 

1.48x10’ 

0.0 

lower 


In Curry et al. 21 , a good experimental correlation was found between the ratio of the 
longitudinal stiffness of the thick to the thin sections of a dropped-ply laminate and the 
compressive strength ratio of the thin section to the dropped-ply laminate (tested 
independently). Despite the fact that this correlation neglects the change in failure mode 
from that of the thin section (compressive strength failure) to that of the dropped-ply 
laminate (delamination), the agreement was very good. A similar approach is now taken in 
looking for correlation between the magnitude of the stiffness change and the interlaminar 
stresses introduced. 

In Cunry’s examination of the effect of stiffness change, strain gage measurements 
were used to determine the average stiffness of the thick and thin sections of the laminate 
and a ratio of the two values was computed. The present study, being purely analytical, 
will rely on an estimation of this ratio based on classical lamination theory (CLT), which 
has already been shown to be in close agreement with the homogenization scheme 
implemented in these analyses for flat laminates. Since the loading in the present study 
consists of an applied N x with e, constrained to zero through the prescription of the out of 
plane deformations, the stiffness coefficient from CLT of interest is A,,. That is, this 
loading can be thought of as an applied strain field of e x prescribed with £,. = 0. Because 
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both symmetric angle-ply and quasi-isotropic laminates have zero shear-extension coupling 
terms A 16 and A 26 , no is induced and the applied N x is related to £° through the relation 
N x =A n e° x in CLT. 

Rather than correlating the stiffness change to a measured failure load as did Curry et 
al., the present study will use the peak value of the delamination fraction. The values of the 
longitudinal stiffness ratio for the five dropped-ply laminates examined are presented in 
Table 6.4, a plot of a linear least squares fit of these ratios vs. the peak 4f values 
presented in Table 6.3 are shown in Fig. 6.12. Although not compared to a baseline case of 
a uniform specimen as did Curry et al., since this would correspond to a zero VF , these 
results indicate that if failure is based on the magnitude of the peak VF , an increase in the 
longitudinal stiffness ratio will result in a proportional decrease in failure load of the 
dropped-ply laminate. 

Table 6.4 Ratios of the longitudinal stiffness of the thick section to the thin section for 
different values of 6 of the dropped sublaminate [±6\. 


e 

(Ai )*«*/( A i),^, 

0° 

1.571 

30° 

1.350 

45° 

1.190 

60° 

1.090 

90° 

1.051 


A similar examination of the in-plane shear loaded cases can be made by using the A ^ 
stiffness coefficient instead of A n . This choice is easily justified by the aforementioned 
decoupling of the shear and extension for symmetric angle-ply and quasi-isotropic 
laminates. That is, since A 16 and are zero, an applied N v induces no £° or £° in CLT, 
which is consistent with the out of plane normal deformation and -^-direction boundary 
condition of the present analysis. Therefore the applied N v is related to through the 
relation N xy =A 66 y° xy in CLT. 
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peak F l/2 



Fig. 6.12 Plot of the correlation between the stiffness ratio and the peak value of V"F 

Gongitudinal compression loading). 


The values of the shear stiffness ratio for the five dropped-ply laminates examined are 
presented in Table 6.5, a plot of a linear least squares fit of these ratios vs. the peak Vf 
values presented in Table 6.3 are shown in Fig. 6.13. 

Table 6.5 Ratios of the shear stiffness of the thick section to the thin section for different 
values of 0of the dropped sublaminate [±0] s . 


e 


0° 

1.076 

30° 

1.337 

45° 

1.425 

60° 

1.337 

90° 

1.076 
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Fig. 6.13 Plot of the correlation between the stiffness ratio and the peak value of VF 

(in-plane shear loading). 


While the correlation here is not as good as for the longitudinal compression case, the 
trend is clearly represented. This is a somewhat curious result because unlike the 
longitudinal compression case, which has significant contribution from both o nn and <T m , 
the value of VF for the in-plane shear loading is entirely due to one stress component, 

The most obvious deficiency here lies in the fact that is symmetric with respect to 6 
about 45° for these layups while the peak values of Vf arc not. The cause of this lack of 
symmetry lies in the difference in transverse shear stiffnesses G 13 and G 1^ for a 
unidirectional ply, which do not enter into CLT. Because of the generalized plane 
deformation assumption of the present theory, the shear strain produced by the applied 
shear loading is due entirely to Bvj dx as Buj By is zero. As a consequence, the load is 
transferred to the dropped plies only through as mentioned above. Therefore, a ply 
oriented such that it has a greater out-of-plane shear stiffness in the y-z plane will have 
larger interlaminar stresses for this shear loading. This would correspond to a 6 closer to 
90° rather than 0° and indeed these are the layups with the larger value of Vf in Fig. 6. 13. 

This difference in G 13 and G 23 also affects the load transfer for the longitudinal 
compression loading, however, in that case both o„ n and G m are major contributors to the 
peak value of the delamination fraction. Thus, the influence of the transverse shear 
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stiffnesses Gu and G 2 3 on the delamination fraction in longitudinal compression is not as 
important as it is in the shear loading case. 

A few other points regarding these correlations are worth noting. First, the fact that 
neither linear least squares fit passes through VF = 0 for a stiffness ratio of 1.0 indicates 
that other influences are contributing to the value of Vf. These other influences are the 
geometry of the taper and the presence of the material discontinuity, neither of which need 
be accompanied by a stiffness change. For example, a laminate with the same taper 
geometry could be made with the dropped plies replaced by resin, leaving a laminate with a 
ratio of CLT stiffnesses of approximately one. A material discontinuity also need not be 
accompanied by a stiffness change as defined here. Wisnom’s' 2 -^ work examining 
laminates with cut, but not discontinued, internal plies is an example of such a 
configuration. As stated in Sec. 1 . 1 . 1 , he found this to be the greatest single factor 
contributing to delamination of unidirectional laminates containing dropped plies loaded 
longitudinally. 

For the longitudinal compression loading, extrapolation of the linear curve fit to a value 
of one for the stiffness ratio would give a delamination fraction of about 1.4 x 10~ 3 while 
the range of values of Vf from (Ai), Aict /(Ai) rtol = 10 to its maximum of 1.57 is only 
about 0.8 x 10 3 . This indicates that these other influences are predominant in contributing 
to VF for this load case. The shear loaded case, on the other hand, would have 
Vf = 1.2 x 10 -3 at a stiffness ratio of one while its range is approximately 2.5 x 10“\ 
indicating that the stiffness change is the major influence. This difference between the 
longitudinal and shear loaded cases also likely contributes to the difference in the quality of 
the correlation. 

In any event, while better in the longitudinal loading case, both of these correlations can 
be useful in rough estimates of a dropped-ply’s likelihood of delamination. More accurate 
predictions, however, require more detailed analyses such as the ones developed in the 
present work. 

6.3 Effect of a Soft Insert 

This parametric study examines the effect of the inclusion of a soft insert ahead of stiff 
dropped plies. The baseline laminate chosen for this study is the same as the asymmetric 
cases with four 0° plies dropped analyzed in the previous two parametric studies and also in 
Sec. 5 . 5.2 where a comparison was made to the finite element results of Curry et al. 21 This 
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laminate is composed of a 20-ply [(±45 / 0 / 90) s / 0 2 ] s from which the center four 0° plies 
are dropped. The insert cases that are compared to the baseline contain four 90° plies 
between the four 0° plies and the thickness change (see Fig. 6.14). The dimension a in 
the figure is the distance from the original location of the end of the 0° plies (x = 0) to the 
new location. Manufacturing such laminates has shown a tendency for the 90° fibers to 
migrate into the area previously modeled as a resin region. Thus, the region which was 
occupied by the neat resin is assumed to be filled with 90 material. 


[±45/0/90] , 


upper interface 


[04]t 


[90 4 ] T 


[±45/0/90] 


Z 


lower interface 


Fig. 6.14 Schematic of dropped ply laminate with soft insert. 

The model used to analyze this problem has the same distribution of six mathematical 
layers as in the previous two studies and in Sec. 5.5.2. The boundary conditions employed 
are the same as those in Sec. 5.5.2 which were dictated by the test conditions used by 
Curry et al. The loading is uniaxial compression with the average normal applied stress in 
the thin section of the laminate designated G m . The dimensions are again normalized by the 
thickness of the dropped plies, t d , which is 0.5588 mm for this laminate. 

The experimental results obtained by Curry et al. for the baseline case show a strong 
tendency for delamination originating at x = 0 either along the top or the bottom of the 
dropped 0° plies. Therefore, the results presented here will be for the interlaminar stresses 
along these two contours, designated the upper and lower interfaces respectively. Figs. 
6.15-6.18 are plots of the interlaminar normal and shear stresses in the contour following 
reference t-n for the baseline case and two soft insert cases. The first soft insert case, 
designated a = 0 (90s) , simply replaces the triangular resin region of the baseline case with 
90° material. The second soft insert case plotted is for a = 4.5^ which is the longest insert 
analyzed. On a side note, the models with the soft insert produce a junction between the 
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Fig. 6.17 Interlaminar normal stresses along lower interface. 
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Fig. 6.18 Interlaminar shear stresses along lower interface. 

From these plots it is clear that inclusion of a soft insert significantly reduces the peak 
interlaminar stresses for this particular laminate. This effect is consistent with what would 
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be expected when the magnitude of the stiffness change at the end of the four 0° plies is 
reduced. Fig. 6.19 shows the effect of the length of the 90° insert on the peak values of the 
tensile interlaminar normal stress and the absolute value of the shear stress. For the baseline 
case, designated w/o insert, the peak tensile interlaminar normal stress occurs at jc = 0 on 
the upper interface and the largest interlaminar shear stress occurs at x = 0 on the lower 
interface. 



a/t d 

Fig. 6.19 Peak interlaminar stresses vs. length of 90° insert. 


It is evident in Fig. 6.19 that the soft insert does not have to be very long to 
significantly reduce the likelihood of delamination initiating at this ply drop-off. Beyond 
filling the resin region with 90° material, which reduces the tensile interlaminar normal 
stress by 39%, the shear by 15%, and the delamination fraction Vf by 29%, little 
improvement is achieved by increasing the length of the soft insert. In fact, in going from 
a = 0 to a = t d the normal stress actually increases slightly. 
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Chapter 7: Concluding Remarks 


7.1 Discussion and Conclusions 

Summary of Approach 

The objective of this research was to develop a stress-based method of approximation 
for the prediction of interlaminar stresses in the vicinity of a ply termination for a composite 
laminate. The aim was to accurately model the geometric and material properties of the 
dropped ply laminate such that interlaminar stresses predicted from the analysis could be 
employed with confidence in a delamination initiation criterion. 

Based on Pagano’s 45 48 - 49 approach for the analysis of axisymmetric involute bodies 
of revolution, the approach chosen modeled the laminate by a series of layers, which can 
have curved boundaries, with the stress field assumed within each layer. Developed for the 
two-dimensional problems of generalized plane deformation and axisymmetry, the stresses 
were assumed to be explicit functions of the thickness coordinate with coefficients , or 
stress variables, that were functions only of the longitudinal coordinate x. The dependence 
of the stresses on the thickness coordinate was determined from the three-dimensional 
elasticity equations of equilibrium. 

The assumed stress field was then substituted into the Hellinger-Reissner variational 
principle, which, after definition of weighted displacement variables in the thickness 
coordinate, was integrated through the thickness. Traction and displacement continuity 
conditions between layers were explicitly imposed on the dependent variables and their 
variations. Stationarity of the functional with respect to all admissible stress and 


132 


Concluding Remarks 



displacement variations consistent with these interlayer continuity conditions leads to the 
governing, or Euler, equations and boundary conditions for the model. Satisfaction of 
these governing equations by the stress and displacement variables, therefore, implies that 
the stress field satisfies equilibrium as well as traction reciprocity between layers and that 
the displacements satisfy continuity conditions between layers. The 
compatibility/constitutive conditions are satisfied in an integral sense within the layers. 

Improvement in Solution Strategy 

The governing equations resulting from the principle constitute a system of both 
differential and algebraic equations, or DAEs. Using Pagano’s 48 - 49 two-step finite 
difference method of solution to this system of DAEs resulted in a solution exhibiting 
unacceptable oscillations when mixed boundary conditions were prescribed. (Pagano only 
solved prescribed traction boundary value problems in Refs. 45, 48, and 49). In the two- 
step finite difference solution, Pagano needed special treatment of the boundary points in 
developing what he called “end” conditions. This special treatment was based on assuming 
the integrand of the functional vanished at the end points, as well as in the open interval 
between the end points. A more detailed evaluation of the system of DAEs allowed a 
resolution of apparent inconsistencies in the number of differential equations and boundary 
conditions without resorting to the use of “end” conditions. This was done through a 
manipulation of some of the equilibrium equations and traction prescribed boundary 
conditions, and resulted in a system with an equal number of differential equations and 
boundary conditions (see Sec. 3.6). The combination of these manipulations and a switch 
to one-step finite differences was found to resolve the deficiency in the numerical solution 
following Pagano’s approach. In addition to the elimination of the oscillations, these 
changes from Pagano’s approach produce a more consistent system mathematically because 
the boundary conditions are distinctly separate from the field differential equations. This 
allows for their solution by other means besides one-step finite differences if desired, while 
Pagano’s interpretation of “end” conditions is specifically tailored to solution by two-step 
finite differences. 

Alleviation of Numerical Ill-Conditioning 

Attempts to apply this new solution strategy to laminated shells of revolution revealed 
severe matrix conditioning difficulties for thin layers at even moderate radii. Development 
of a second model based on the assumption of generalized plane deformation in cartesian 
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coordinates was undertaken to determine the cause of the conditioning. Initially the model 
was developed using the same approach to deriving the governing equations as Pagano 
used, and conditioning problems were present for layers located at moderate values of the 
thickness coordinate. While providing some benefit, a series of changes to the stress 
assumptions including scaling of the bubble functions and imposing orthogonality on them 
were unsuccessful in overcoming the problem. The change that had the greatest positive 
impact on the conditioning difficulty was the modification of the displacement weighting 
functions from simple powers in the thickness coordinate to functions based on the shape 
functions used for the stresses. These modifications to the stress assumptions and weighted 

displacements are detailed in Secs. 3.2 and 3.3. 

These changes allowed the generalized plane deformation model to handle problems 
with very thin layers, as would be encountered in modeling of laminates with dropped 
plies. However, when these remedies were applied to the axisymmetric model, the results 
were not as good. While the conditioning was greatly improved, patch test problems 
revealed that numerical conditioning was marginal even for layers that were thick relative to 
those that would need to be modeled in the analysis of cylindrical shells with axially 

dropped plies. 

Sensitivity to Solution Accuracy 

Attempts at solutions of these models by trapezoidal one-step finite differences revealed 
that they are very sensitive to the accuracy of the solution strategy employed. Unable to 
produce acceptable results for problems with rigid body displacements, the 0(h 2 ) accurate 
trapezoidal scheme was replaced with the OQt) accurate two-stage Gauss implicit Runge- 
Kutta scheme. This improvement in accuracy did not come without cost, however. 
Implementation of the Gauss scheme required conversion of the algebraic equations to 
differential and the inversion of a 2n x 2n matrix, n being the order of the system, at each 
finite difference step. In addition, this approach to solution was found to cause a 
degradation in solution matrix conditioning as compared to the application of the trapezoidal 
scheme to the DAEs. This degradation in conditioning was found to be caused mostly by 
the conversion of the DAEs to a fully differential system. 

Reformulation of Axisymmetric Model 

Along with the changes made to the axisymmetric model for the purpose of improving 
the reliability of the solutions obtained, a further modification was made to Pagano s 
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original approach. Through careful selection of the dependence of the “in-plane” stress 
components on the radial coordinate, the number of stress variables required to allow 
satisfaction of equilibrium was reduced by one per layer as was the number of displacement 
variables per layer. While this is a modest reduction in the number of unknowns per layer 
(from 25 to 23), it is significant and actually required less effort in the development of the 
model as compared to Pagano’s. It also provided the added benefit of having a one-to-one 
correspondence in the stress and displacement unknowns to the generalized plane 
deformation model, allowing a simpler conversion in the coding of the two models. 
Unfortunately, the difficulties in obtaining accurate and reliable solutions with the 
axisymmetric model for problems subjected to internal pressure did not allow these 
advances to be utilized. 

Modeling Multiple Plies in a Layer 

A homogenization scheme was applied in order to allow multiple plies to be represented 
by a single layer in one of the models. This scheme used a combination of the Voigt and 
Reuss approximations in deriving compliance coefficients for a single homogeneous layer 
equivalent in a strain energy sense to the sublaminate comprising the actual layer. This 
approach to modeling multiple plies in a single layer was compared to piecewise evaluation 
of the weighted integrals of the compliances and was found to be superior due to a lack of 
in-plane displacement continuity enforcement in the latter approach (see Sec. 5.2.3). 

Verification Studies 

Patch test problems were used to evaluate the accuracy and limitations of the models 
and their solution by the finite difference approximation. The term patch test is used 
because of the similarity to tests of the same name used to verify the accuracy of finite 
elements. Essentially, a problem is chosen to which the exact elasticity solution is known, 
which is usually a state of uniform stress and strain. By geometrically distorting the layers 
subdividing the body, some limits on the accuracy of the numerical solution method can be 
determined. These patch tests revealed the difficulties with the axisymmetric model 
mentioned previously. They also helped establish the level of confidence in the generalized 
plane deformation model necessary to move on to the analysis of laminated composites. 

Evaluation of the accuracy of the generalized plane deformation model for laminated 
composites was also done through comparing solutions obtained from it to established 
solutions. The first such comparison was to the tensile coupon free edge problem which 
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Pagano 44 compared to the highly refined finite element solution of Wang and Crossman. 57 
The results for this case showed near perfect agreement for the interlaminar stresses. While 
this problem did not involve the geometric complexity of a dropped ply, the results 
indicated that the model was capable handling the stress gradients present in a region of 
stress concentration. 

The second solution compared to was that of Curry et al. 71 and their finite element 
analysis of a compression loaded asymmetric dropped-ply test specimen. This extended the 
evaluation of the model to an examination of it’s ability to model complex geometries with 
very thin layers as well as stress concentrations. While the interlaminar stresses predicted 
by the two approaches did not match as well as the free-edge problem, the correlation was 
good. This showed that the model had the ability to model complex problems such as 
dropped-ply laminates. 

Effect of Dropped-Ply Laminate Eccentricity 

Parametric studies were done to evaluate some of the factors affecting the stress 
concentration in dropped-ply laminates. Eccentricity of the middle surface was examined 
for a dropped-ply laminate composed of continuous quasi-isotropic sublaminates of 
[±45 / 0 / 90] t layup, and with dropped sublaminates of either [0 4 ] x , [90 4 ] T , or [±45] s 
layups. Material properties for AS4/3502 graphite/epoxy were used. Two geometries were 
evaluated for each of the three layups of the dropped sublaminates, one with eccentricity 
and one without Each of these configurations was loaded in either longitudinal 
compression, transverse tension, or in-plane shear. The relative propensity of each case to 
delaminate was evaluated based on a quadratic delamination criterion using the interlaminar 
stresses. For the longitudinal compression and in-plane shear load cases, eccentricity was 
determined to have a minimal effect on the likelihood of delamination occurring. The 
transverse tension case was not considered likely to experience delamination failure at all 

due to the ply drop-off. 

Effect of Stiffness Discontinuity 

The eccentricity study clearly indicated that the stiffness discontinuity introduced by 
dropping plies had a much greater influence on the interlaminar stresses than did 
eccentricity. A second parametric study was performed to examine the correlation between 
the magnitude of the stiffness discontinuity and the likelihood of delamination for 
asymmetric dropped-ply laminates. The stiffness of the dropped sublaminate was varied 
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through changing the ply angle 6 in a [±0] s sublaminate layup. The material used was 
again AS4/3502 graphite/epoxy, and the loadings examined were longitudinal compression 
and in-plane shear. For the longitudinal compression loading, a good correlation was found 
between the delamination criterion’s prediction of the likelihood of delamination and the 
ratio of the longitudinal stiffness of the thick section of the laminate to that of the thin 
section. A similar correlation for the in-plane shear loaded cases using the shear stiffness 
discontinuity instead of the longitudinal stiffness was not quite as good, but did 
approximate the trend acceptably. 

Influence of Soft Insert 

A third study of dropped-ply laminates examined the effect of inserting “softer” 90° 
material ahead of four 0° dropped plies in an effort to reduce the peak interlaminar stresses 
under longitudinal compression. AS4/3502 graphite/epoxy material was used along with 
the same asymmetric geometry of taper as used in the eccentric case. The length of the soft 
insert was increased incrementally to determine an optimal length. It was found that simply 
filling the region ahead of the dropped plies previously occupied by resin had the most 
significant effect on the interlaminar stresses, reducing the interlaminar normal by nearly 
40%. Beyond that, increasing the length of the insert provided only marginal benefit. 

7.2 Suggestions for Future Work 

At present, a major weakness of these models and their solution lies in the large 
computational burden in both storage and CPU time. This burden stems from three areas; 
the use of a banded matrix solver for a matrix with block structure, the large number of 
dependent variables, and the high degree of accuracy required of the solution method. 
Compounding this problem is the issue of conditioning of the solution matrix, which 
degraded with the application of the more accurate Gauss implicit Runge-Kutta scheme due 
mainly to the conversion of the DAEs to a fully differential system. Therefore, 
recommendations for future work on these models should begin with improvements in their 
solution with the easiest being the use of a block diagonal solver. This was not done in the 
present work because such a solver that also estimates condition numbers was not 
available. Elimination of variables that appear only algebraically in the governing equations 
through manipulation could reduce the number of unknowns per layer in the finite 
difference solution by nine. However, this procedure could likely only be done numerically 
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in a practical sense, and therefore would involve matrix inversions. The cost of these 
inversions, which would reduce the computational benefits of this approach, must be 
considered as well. Another approach worth consideration is adapting a higher order 
accurate one-step finite difference scheme, such as the Gauss scheme used here, to the 
solution of the original DAE form of the governing equations. This would avoid the 
conversion to a fully differential system, which was found to degrade the numerical 
solution conditioning significantly. Finally, some form of finite element solution such as 
the one developed by Sandhu et al. 50 deserves more attention than was possible in the 
present study. 

These improved solution strategies should be implemented first for the generalized 
plane deformation model. Those that are most successful could then be applied to the 
axisymmetric model in an attempt to anive at accurate solutions for laminated cylindrical 
shells subjected to internal pressure. 

Other improvements to the models worth considering include the development of a 
more refined homogenization scheme (see Sec. 5.2 and Appendix A) that would provide 
compliance coefficients that vary through the thickness of a layer. This would account for 
some of the effect of the stacking sequence within a layer on the compliances and may 
produce more accurate results with fewer layers. Finally, developing continuity conditions 
that would allow two layers to be joined to form a single layer would provide the dual 
benefit of reducing the number of layers where they are not needed and minimizing the 
number of very thin layers in a model of such features as ply drop-offs. 
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Appendix A: Homogenization 
Procedure 


This appendix presents the homogenization procedure applied to the mathematically 
modeled layers that contain multiple plies (a sub-laminate). It is based on the equivalence of 
the strain energy between a representative volume of the heterogeneous sub-laminate, for 
which the material properties are known on a constituent basis, and a representative volume 
of homogeneous solid, for which the material properties are desired. In equating the strain 
energies, a combination of the Voigt 58 and Reuss 59 approximations are applied such that 
in-plane compatibility and out-of-plane equilibrium are enforced between plies. 

A.l Equating Strain Energy 

A schematic of representative volumes of the heterogeneous and the homogeneous 
materials is shown in Fig. A. 1 . Both elements have in-plane dimensions Ax and Ay and 
total thickness t. The heterogeneous element may be composed of any number of plies, n, 
with each ply having arbitrary thickness f, within the restriction that = t. Therefore, the 
volume fraction of each constituent ply is v, =tjt. In-plane dimensions Ax and Ay are 
assumed to be infinitesimal, and the stresses and strains within a ply are assumed spatially 
uniform. Each ply is assumed linear elastic, homogeneous, and anisotropic. 
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Fig. A.l Heterogeneous sub-laminate and equivalent homogeneous solid. 

The linear-elastic strain energy for the two volume elements are, for the heterogeneous 
sub-laminate: 


(A-l) 

^ i=l 

and for the equivalent homogeneous solid: 

^0 k e k tAxAy (A-2) 

where o k and £* are the stresses and engineering strains in contracted notation, i.e. 
[o x ,a y ,a z ,<T yz ,a xzJ o xy ] T = [ <r, , <r 2 , , cr 4 , cr 5 , a 6 ] r with the same order for the engineering 
strains. Repeated subscripts on stress and strain components are summed out from one to 
six in the usual indicial notation. In Eq. (A-l), the superscript in parentheses denotes the 
layer that the stress or strain component is associated with. The stresses and strains for the 
homogeneous solid have no superscripts. Equating these expressions, eliminating like 
terms and employing the definition of the volume fraction leaves 

( A ‘3) 

i=i 

Substituting Hooke’s Law, expressed as 

o k = C w £, or e k = S kl o l (A-4) 

where C and S are the stiffness and compliance matrices respectively, the equivalence of 
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strain energy, Eq. (A-3), can be written in the following forms: 


C w £*£ ; = 

i = l 

(A- 5a) 

s u o k o, = 2v^v; ) o{ , ' > 

(A-5b) 


i=l 


in which C® and S® denote the known stiffness and compliance matrices for the ith ply. 
The goal of the remaining derivation, therefore, is to determine the C and S matrices 
occurring on the left sides of Eqs. (A-5), i.e. the homogenized stiffness and compliance 

matrices. 

A stress concentration matrix K (l) for the ith ply, which relates the stress in the 
homogenized composite to the stresses in the individual constituent plies, is defined by: 


— (0 _ 

°/fc ~ A W °l 


(A-6) 


Substituting into Eq. (A-5b) yields 




i=l 


for every state of stress in the homogeneous solid. Therefore, 


(A-7) 


1=1 


(A-8) 


and determining the homogenized compliance matrix S reduces to the determination of the 
stress concentration matrices K® and the application of Eq. (A-8). 

Similarly, a strain concentration matrix L® is defined by 


,(i) _ r( 0 


= 


(A-9) 


Substituting into Eq. (A-5a) leads to 


n 



i=l 


(A- 10) 


So determination of C reduces to the determination of the strain concentration matrices L®. 
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A.2 Voigt and Reuss Approximations 


The stress concentration matrix is determined through application of a combination of 
the Voigt 58 and Reuss 59 approximations. These approximations were developed for the 
micromechanical evaluation of average elastic moduli of materials with multiple constituents 
such as polycrystals. The Voigt approximation assumes that the strain in each constituent is 
uniform and equal to an average strain. The Reuss approximation assumes that the stress in 
each constituent is uniform and equal to an average stress. It has been proven by Hill 60 that 
the Voigt approximation gives upper bounds and the Reuss approximation gives lower 
bounds for determination of effective elastic moduli for a composite. 

The Voigt approximation is applied to the in-plane strains and serves to enforce in-plane 
displacement compatibility between plies. The Reuss approximation is applied to the out- 
of-plane stresses and serves to enforce out-of-plane equilibrium between plies. These 
conditions are stated as 


£ — 
c, c, 

e 2 =£ 2 ) I = 1,2,.. .,71 

t 6 


and 


(A-l 1) 



c 5 = of 


(A- 12) 


In addition, conditions are placed on the out-of-plane strains and the in-plane stresses. 
The conditions on the out-of-plane strains are stated as 


t£i = X'.ef 

te 4 = J>f < A ' 13 ) 

re 5 = £f,ef 

and are imposed in order to equate the relative displacements between the lower and upper 
surfaces of the volume elements due to the out-of-plane strains. The conditions on the in- 
plane stresses are stated as 
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(A- 14) 


= i (i) 

^2 = X^'’ 

^6 = 

and are imposed in order to equate the in-plane forces on the vertical faces between the 
volume elements. Division through by the thickness t gives the form for Eq. (A- 13) of 

= X V ' e 3 0 

£4 = X V > £ 4 ’ (A ‘ 15) 

£5 = X V * e 5 ° 

and for Eq. (A- 14) 

<*1 =X v ^i <0 

cj 2 = Xv,^'* (A- 16) 

^6=X V - <T 6 ) 


A.3 Homogenized Compliances 


Based on Eqs. (A-6), (A- 12), and (A- 16), the structure of the stress concentration 
matrix for layer i has the form 


*n 

*,2 

*13 

*14 

*15 

*16 

*21 

* 22 

*23 

*24 

*25 

*26 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

*6. 

*62 

*63 

*64 

*65 

*66 


(A- 17) 


Determination of the 18 terms that are unknown is accomplished through the application of 
the conditions of Eq. (A- 11) along with Hooke’s law, the in-plane equilibrium conditions 
of Eq. (A- 16) and the out-of-plane equilibrium conditions of Eq. (A- 12). 

Because the three in-plane strains for the homogeneous composite do not enter into the 
determination of the stress concentration matrix, they are eliminated from Eq. (A-l 1) 
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leaving 3(n-l) conditions: 


r (0 - pO+n 
c i "" c \ 

£ 2 ’ = e 2 +1> / = l,2,...,n- 1 (A- 18) 

4° = ^‘ +l) 

Substitution of the strain-stress form of Hooke’s law into Eq. (A- 18) and accounting for 
the out-of-plane equilibrium conditions of Eq. (A- 12) yields 

SfJ'of ‘ + Si‘W 6 ° - Sj{ +1) of i+1) - Otf +1) “ S iT I>cr 6 +1> 

= (sir 1 - + w - sfi’K 


c(i )—.(») , nO'),— (i) , r(* )—.(«) __ ofi + D — .(i+U _ r(i+l) — »+l) _ r(*+l) —.( j+1) 
°21 °1 ^ ^22 °2 °26 °6 °21 °22 U 2 °26 U 6 

= ($£“’ - + (s;;* 11 - s£)o, + (s"* 11 - ^>)a 5 


3X’ + ago?’ + s“>o"’ - sjrx* 11 -s^V“> - 
= (C - «>, + (&“’ - S£’K + (sT 1 -ag’)®, 

for i = l,2,...,n-l. 

Defining the 3 x 3 matrix of in-plane compliances for the ith ply as 

Sji s 12 s 16 

S 2[ S 22 S 26 

$bl $62 ^66 

as well as the 3 x 3 matrix of differences in the compliances that couple in-plane and out- 
of-plane response as 

" ofi+l) r>(i) r(i+l) _ c(0 cO’+l) c(0 

°13 °13 °14 °14 °15 °\5 

o(0 — cd+U c(0 c<*+i) c(0 cO'+l) __ o(0 /a 

O — 0 23 0 23 0 24 024 025 0 25 v™ 

r(i+t) _ r(i) r(i+ 1) _ c(0 c(i+l) _ r<0 

_°63 °63 °64 °64 °65 °65 _ 

and the 3 x 1 vector of in-plane stresses as 



<o 

(A-20) 
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<T (,) = < 


(i) 


(A-22) 


Eqs. (A- 19) can be rewritten in matrix form as 

S <0 a <i) -S <i+1) CT (i+u = S (,) <{ 


(A-23) 


which again is for i = 1, 2, ... , n-1. By defining the 3 x 3 volume fraction matrix as 


V (i) = 


v, 0 0 

0 v, 0 

0 0 v, 


(A-24) 


Eq. (A- 16) can be rewritten in the form 

V (1) {f <l) + V (2) cf (2) + • • ■ + V <B) CT <n) = ^ 


l a e) 


(A-25) 


Eqs. (A-23) and (A-25) can be combined into one matrix equation: 


r^on 

<t (2) 


r(n) 


= B 


1°5J 


(A-26) 


The partitioned form of A, a 3 nx 3 n matrix, is 
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' yd) 

y<2> 

y< 3 > 


y(n) 

s (,) 

_S<2) 

0 

0 

0 

0 

s <2) 

-s <3) 

0 

0 

0 

0 

s <3) 

0 

0 

0 

0 

0 

s ( "‘ l) 

1 

tzy .. 

i 


and the partitioned form of the 3n x 6 matrix B is 


B = 


I 

0 

0 


0 

s (1) 

s (2) 


0 s (n_1) 


(A-27) 


(A-28) 


In order to simplify the partitioned nature of the A and B matrices, the vector of 
homogenized stresses on the right hand side of Eq. (A-26) was reordered. In defining the 
stress concentration matrices K (,) in Eq. (A- 17), the conventional order was used. 
Therefore, a transformation matrix T is defined as 


'1 0 0 0 0 0 

0 1 0 0 0 0 

0 0 0 0 0 1 

0 0 1 0 0 0 

0 0 0 1 0 0 

0 0 0 0 1 0 


(A-29) 


and solving for the heterogeneous in-plane stresses in Eq. (A-26) yields 


r^on 

a {1) 


a in) 


> = A -1 BT< 


l <*6 


(A-30) 


Therefore, the undetermined portion of the stress concentration matrices, defined by 
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K ,0 = 


*n 

K 

K, 


21 


61 


are determined by 


K n K n 

^22 ^23 

*62 *63 


K (1) 

K (2) 


*,4 *» 

* 24 * 25 
k m * 65 


* 


i(.) 


16 


K. 


26 


K, 


66 


(A-31) 


K 


(*) 


= A BT 


(A-32) 


and the homogenized compliance matrix is defined by Eq. (A-8). 

A.4 Homogenized Stiffnesses 

Eqs. (A-9), (A- 11), and (A-15) imply the structure of the strain concentration matrix in 
the ith ply is 

(i) 


L <0 = 


Eqs. (A- 12) imply 


1 

0 

0 

0 

0 

0 ' 

0 

i 

0 

0 

0 

0 

^31 

Lj 2 

^33 

*34 

*35 

636 

4l 

L 42 

^43 

^44 

*45 

*46 

*3, 

^52 

^53 



*56 

0 

0 

0 

0 

0 

1 

4° = 

<Tj +1) 





W = 

of 1 ’ 

i 

= 1,2,. 

..,n- 

1 


(A-33) 


(A- 34) 


of = cr^ 0 


Substitution of the stress-strain form of Hooke’s law into Eq. (A-34) and using Eqs. (A- 
1 1) leads to 
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£(i)—{i) _ £(i + l)£(i+l) _ £»(i) 


L £ 6J 


i = 1,2 - 1 


(A-35) 


in which 


C (i) = 


Qj ^35 


c c c 


45 


43 ^44 

53 £54 Qs . 


(i) 


(A-36) 


e (,) = 


l £ 5 


(i) 


(A-37) 


and 


C (i> = 


c5i*”-c: 


"41 
'M i + 

"51 


WO 

Hi 

C 32 

WO 

H 2 

Wi+1) 

^36 

WO' 

C 36 

wo 

Hi 

Wi+1) 

H2 

WO 

H2 

Wi+1) 

H6 

WO 

°46 

wo 

Hi 

Wi+1) 

H 2 

WO 

H 2 

/■^(i+l) 

He 

Wi) 

H6 _ 


(A-38) 


In matrix form Eq. (A- 15) is 


yU)£<i> + v (2) e <2) + • ■ • + V tB, £ ( "’ = 


(")?(«) _ 


'5J 


(A-39) 


where matrices V (i) are given by Eq. (A- 24). Eqs. (A-35) and (A-39) form a 3n x 3n 
system of equations for £ <0 given by 


A' 


£ (2) 


£ <b) 


= B' 


l £ 6 


(A-40) 


in which matrix A is the same form as A in Eq. (A-27) with submatrices S' 1 replaced by 
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C (,) , and matrix W is the same form as B in Eq. (A-28) with S <n replaced by C <0 . Let the 
nontrivial three rows of L (i) be denoted by 


£31 ^32 ^33 ^34 ^35 ^36 

L 4] L 42 L 43 L 45 

^51 ^52 ^53 ^54 ^55 ^ 56 . 


Then it can be shown that 


'C ir 

l (2) 


(AT'B'r 


L <») 


where 


'0 

0 

1 

0 

0 

O' 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 


(A-41) 


(A-42) 


(A-43) 


Recently, Sun and Li 61 have presented explicit formulas for the homogenized elastic 
moduli of a sub-laminate. Their formulation also uses Eqs. (A-ll), (A- 12), (A- 15), (A- 16) 
and Hooke’s law as presented here, but they do not use the stress and strain concentration 
matrices, Eqs. (A-6) and (A-9), and equivalence of the strain energy as represented by Eq. 
(A-3). 
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Appendix B: Two-Stage Gauss 
Implicit Runge-Kutta Finite 
Difference Scheme 


This appendix presents the two-stage Gauss implicit Runge-Kutta finite difference 
scheme. This is a finite difference adaptation of a subclass of the Runge-Kutta schemes for 
implicit integration of differential equations as presented on pp. 210-216 of Ref. 52 and 
summarized here. It should be noted that the notation used here is the same as that of Ref. 
52 and that indicies or variables that coincide with previous chapters do not imply any 
relationship between the two. 

B.l Implicit Runge-Kutta Schemes 

The Jfc-stage Runge-Kutta scheme for the system y' = f(x,y) over N intervals is given 
by 

y i+ i = y i +h i^Pj f ij \<i< n (B-i) 

j=i 

where 


< B - 2 > 

/=i 

are the evaluations of the function f at the points 
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Xij =x,+ hfij 


1 <j<k, 1 <i<N 


(B-3) 


which are located by the values of p jt with 

0 <a <p 2 <-<p t <1 (B- 4 ) 

Subscript i denotes a mesh interval, while subscripts; and l are stage counters within a 
particular interval. Points x v are sometimes called collocation points. 

The choice of the values of the parameters p p and Cfy, j,l = defines the 
specific Runge-Kutta scheme. If p, = 0 and a jt = 0 for j < l, then the scheme is explicit, 
while if either condition does not hold, the scheme is implicit. Runge-Kutta schemes are 
typically represented in tableau form as shown in Eq. B-5. 


Pi 

(X u 

- 

Pk 


•" a kk 


A 

... (5 k 


(B-5) 


B.2 Development of Finite Difference Scheme: 
Parameter Condensation 

In order for the Runge-Kutta schemes to be applied in one-step finite difference 
approximations, they must be written in the form 

s,y.+R,y,.,=q, (b-6) 

where y, and y i+1 are the mesh values of what are referred to as global unknowns because 
they are determined through solution over the entire domain. For the Runge-Kutta 
schemes, however, there are also variables local to the interval [x,,x 1+1 ] for each i, 

1 < i < N. The local unknowns for the ith interval are 

< B - 7 > 

Elimination of these local unknowns, referred to as parameter condensation from the 
analogous procedure in finite elements, leaves an expression in the form of Eq. (B-6). 

For the case of a linear boundary value problem, given by 
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y' = f(*,y) = i4(x)y + q(jr) 


(B-8) 


where i4(x) is an nxn matrix, eq. (B-2) becomes 

f a = ^(^)[y, + h^a jt f u h qU,y) 1 (B-9) 

i = i 

Using the notation of Eq. (B-7), Eq. (B-9) can be rearranged and written in matrix form for 
interval i as 


% f, = ^y,+q, 

where 

q(*a) 

|q(-**). 

By defining the matrix D, as 



— <* lk A(x a ) 

V - 

~Mx a y 



* *1 — 



(B-10) 


(B-l 1) 


A/) 

Eq. (B-l) can be restated as 

y, +] = y, + \ d, f, 

which, through substitution of Eq. (B-10) solved for f , can then be expressed as 

y, + i = r ,y,+ r , 

through the definition of T, and r, as 


(B-l 2) 


(B-l 3) 


(B-14) 


(B-l 5) 


Eq. (B-14) is now in the form of Eq. (B-6) with S t = T,, /?, = / and q, = r; . 


B.3 Gauss Schemes 


Gauss schemes are characterized by the choice of points p,, ..., p* as the zeroes of a 
Legendre polynomial. The 2-stage ( k = 2) Gauss scheme used in the present study has the 
tableau 
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(B-16) 


Since p, * 0 , this scheme is implicit and the parameter condensation approach of the 
previous section was used. This particular scheme was chosen because of its implicitness, 
which produces better stability, its high order of accuracy, 0(h 4 ), and its efficiency in 
obtaining this accuracy. 

The major drawback of this approach lies in the potential computational burden of 
inverting the bixkn matrix W, for each interval in the finite difference mesh. This expense 
was reduced in the present work by using a uniform mesh spacing h, as much as possible 
in regions of constant coefficient governing equations, allowing the repeated use of the 
same computed . 
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